Analysis

Author

Dr. Alexander Bauer

Published

August 20, 2024

Code
library(readxl)     # read xlsx data
library(dplyr)      # data handling
library(tidyr)      # data transformation
library(ggplot2)    # data visualization
library(patchwork)  # join ggplot graphics
library(ggpubr)     # ggplot handling
library(kableExtra) # print tables

library(mgcv)     # regression model estimation
library(APCtools) # regression model visualization
library(plotROC)  # ROC curve and AUC calculation

# set global ggplot theme
theme_set(
  theme_minimal() +
    theme(plot.title       = element_text(hjust = 0.5),
          plot.background  = element_rect(fill = "white", color = "white"),
          panel.background = element_rect(fill = "white", color = "white"))
)
Code
# colors for the main areas of analysis
col_vector <- c("knowledge" = "#2F9D64",
                "beliefs"   = "#00868B",
                "barriers"  = "#0062A3")

Data preparation

Code
# read data
dat <- read_excel("../data/Plant-Based Survey - Registered Dietitians_V1.0 .xlsx",
                  sheet = "Combined", skip = 2, n_max = 335)

# read lookup table from the first two rows in the Excel file
dat_lookup <- read_excel("../data/Plant-Based Survey - Registered Dietitians_V1.0 .xlsx",
                         sheet = "Combined", n_max = 2)
Outlier handling

One person (ID 135) has age value 2. Since all other responses by this person seem plausible we assume this was a typo and set the person’s age value to missing.

Code
dat$`What is your age (in years)?`[dat$`What is your age (in years)?` == 2] <- NA
Code
### properly format the lookup dataset
dat_lookup <- dat_lookup %>% 
  t() %>% 
  as.data.frame() %>% 
  rename(group = V1,
         label = V2) %>% 
  mutate(variable = row.names(.)) %>% 
  select(group, variable, label)
row.names(dat_lookup) <- NULL


### format the main dataset
# use the short variable labels for the main dataset
colnames(dat) <- dat_lookup$variable

# properly encode variables
dat <- dat %>% 
  mutate(age_cat              = case_when(age < 25  ~ paste0(min(age, na.rm = TRUE), "-", 24),
                                          age < 30  ~ "25-29",
                                          age < 35  ~ "30-34",
                                          age < 40  ~ "35-39",
                                          age < 45  ~ "40-44",
                                          age < 50  ~ "45-49",
                                          age < 55  ~ "50-54",
                                          age >= 55 ~ paste0("55-", max(age, na.rm = TRUE))),
         age_cat              = factor(age_cat),
         cohort               = factor(cohort),
         survey_consent       = factor(survey_consent),
         registered_dietitian = factor(registered_dietitian),
         gender               = factor(gender, levels = c("Male", "Female")),
         education            = case_when(education == "student"              ~ NA_character_, # unsure about this person
                                          grepl("post", tolower(education))   ~ "Postgraduate degree",
                                          education == "higher diploma"       ~ "Postgraduate degree",
                                          grepl("master", tolower(education)) ~ "Postgraduate degree",
                                          TRUE                                ~ education),
         education            = factor(education, levels = c("Undergraduate degree", "Postgraduate degree", "PhD")),
         dietetics_country    = case_when(tolower(dietetics_country) == "ireland" ~ "Ireland",
                                          TRUE                                    ~ "UK"),
         dietetics_country    = factor(dietetics_country),
         dietitian_forHowLong_cat      = case_when(dietitian_forHowLong < 3   ~ "up to 3 years",
                                                   dietitian_forHowLong < 7   ~ "4 to 6 years",
                                                   dietitian_forHowLong < 10  ~ "7 to 9 years",
                                                   dietitian_forHowLong < 15  ~ "10 to 14 years",
                                                   dietitian_forHowLong < 20  ~ "15 to 19 years",
                                                   dietitian_forHowLong < 25  ~ "20 to 24 years",
                                                   dietitian_forHowLong >= 25 ~ "25+ years"),
         dietitian_forHowLong_cat      = factor(dietitian_forHowLong_cat,
                                                levels = c("up to 3 years",  "4 to 6 years",
                                                           "7 to 9 years",   "10 to 14 years",
                                                           "15 to 19 years", "20 to 24 years",
                                                           "25+ years")),
         dieteticsArea_academia        = grepl("Academia", dietetics_area) %>%         factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsArea_catering        = grepl("Catering", dietetics_area) %>%         factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsArea_community       = grepl("Community", dietetics_area) %>%        factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsArea_hospital        = grepl("Hospital", dietetics_area) %>%         factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsArea_industry        = grepl("Industry", dietetics_area) %>%         factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsArea_primaryCare     = grepl("Primary care", dietetics_area) %>%     factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsArea_privatePractice = grepl("Private practice", dietetics_area) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsArea_publicHealth    = grepl("Public health", dietetics_area) %>%    factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsArea_research        = grepl("Research", dietetics_area) %>%         factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_cardiothoracic   = grepl("Cardiothoracic",          dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_elderly          = grepl("Care for the elderly",    dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_diabetes         = grepl("Diabetes",                dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_gastroenterology = grepl("Gastroenterology",        dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_ICU              = grepl("Intensive care unit",     dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_ltCareFacility   = grepl("Long-term care facility", dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_maternity        = grepl("Maternity",               dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_nutritionSupport = (grepl("nutrition support",       tolower(dietetics_subspecialty)) |
                                                     grepl("enteral",               tolower(dietetics_subspecialty))) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_obesity          = grepl("Obesity",                 dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_oncology         = grepl("Oncology",                dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_paediatrics      = grepl("Paediatrics",              dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_psychEating      = grepl("Psychiatry/eating disorders", dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_renal            = grepl("Renal",             dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_sportsNutrition  = grepl("Sports nutrition",  dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_transplantation  = grepl("Transplantation",   dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_weightManagement = grepl("Weight management", dietetics_subspecialty) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_neurology        = grepl("neuro",             tolower(dietetics_subspecialty)) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         dieteticsSubspecialty_obesityOrWeightManagement = case_when(dieteticsSubspecialty_obesity == "yes"          ~ "yes",
                                                                     dieteticsSubspecialty_weightManagement == "yes" ~ "yes",
                                                                     TRUE                                            ~ "no"),
         dieteticsSubspecialty_obesityOrWeightManagement = factor(dieteticsSubspecialty_obesityOrWeightManagement,
                                                                  levels = c("no", "yes")),
         diet_personal                       = tolower(diet_personal),
         diet_personal                       = case_when(grepl("wfpb", diet_personal)       ~ "WFPB",
                                                         grepl("omnivore", diet_personal)   ~ "Omnivore",
                                                         grepl("flex", diet_personal)       ~ "Flexitarian",
                                                         grepl("pesc", diet_personal)       ~ "Pescetarian",
                                                         grepl("mediterr", diet_personal)   ~ "Mediterranean",
                                                         grepl("vegetarian", diet_personal) ~ "Vegetarian",
                                                         grepl("vegan", diet_personal)      ~ "Vegan",
                                                         grepl("plant", diet_personal)      ~ "Vegan",
                                                         TRUE                               ~ "Omnivore"),
         diet_personal                       = factor(diet_personal, levels = c("Omnivore", "Mediterranean", "Flexitarian", "Pescetarian", "Vegetarian", "Vegan", "WFPB")),
         diet_perceptionPBdiet_DASH          = grepl("DASH diet",  diet_perceptionPBdiet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         diet_perceptionPBdiet_vegan         = grepl("Vegan diet", diet_perceptionPBdiet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         diet_perceptionPBdiet_vegetarian    = grepl("Vegetarian diet",    diet_perceptionPBdiet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         diet_perceptionPBdiet_flexitarian   = grepl("Flexitarian diet",   diet_perceptionPBdiet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         diet_perceptionPBdiet_meditteranean = grepl("Mediterranean diet", diet_perceptionPBdiet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         diet_perceptionPBdiet_MIND          = grepl("MIND diet",      diet_perceptionPBdiet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         diet_perceptionPBdiet_portfolio     = grepl("Portfolio diet", diet_perceptionPBdiet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         diet_perceptionPBdiet_WFPB          = grepl("Whole food plant-based diet", diet_perceptionPBdiet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         diet_perceptionPBdiet_Planet        = grepl("Eat Lancet Planetary Health diet", diet_perceptionPBdiet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         across(starts_with("wpWFPB"), function(x) { factor(x, levels = c("No, strongly disagree", "No, disagree", "Not sure", "Yes, agree", "Yes, strongly agree")) }),
         plantProtein_isIncomplete = factor(plantProtein_isIncomplete, levels = c("No, strongly disagree", "No, disagree", "Not sure", "Yes, agree", "Yes, strongly agree")),
         WFPB_micronutrientsOfConcern_calcium   = grepl("Calcium", WFPB_micronutrientsOfConcern) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         WFPB_micronutrientsOfConcern_choline   = grepl("Choline", WFPB_micronutrientsOfConcern) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         WFPB_micronutrientsOfConcern_folate    = grepl("Folate", WFPB_micronutrientsOfConcern) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         WFPB_micronutrientsOfConcern_iodine    = grepl("Iodine", WFPB_micronutrientsOfConcern) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         WFPB_micronutrientsOfConcern_iron      = grepl("Iron", WFPB_micronutrientsOfConcern) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         WFPB_micronutrientsOfConcern_LComega3  = grepl("Long-chain omega-3", WFPB_micronutrientsOfConcern) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         WFPB_micronutrientsOfConcern_potassium = grepl("Potassium", WFPB_micronutrientsOfConcern) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         WFPB_micronutrientsOfConcern_selenium  = grepl("Selenium", WFPB_micronutrientsOfConcern) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         WFPB_micronutrientsOfConcern_SComega3  = grepl("Short chain omega-3", WFPB_micronutrientsOfConcern) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         WFPB_micronutrientsOfConcern_vitB12    = grepl("Vitamin B12", WFPB_micronutrientsOfConcern) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         WFPB_micronutrientsOfConcern_vitD      = grepl("Vitamin D", WFPB_micronutrientsOfConcern) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         WFPB_micronutrientsOfConcern_thiamine  = grepl("Thiamine", WFPB_micronutrientsOfConcern) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         WFPB_micronutrientsOfConcern_zinc      = grepl("Zinc", WFPB_micronutrientsOfConcern) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_alzheimer        = grepl("Alzheimer", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_heartDisease     = grepl("Heart disease", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_highCholesterol  = grepl("High cholesterol", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_hypertension     = grepl("Hypertension", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_T2Diabetes       = grepl("Type 2 Diabetes", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_obesity          = grepl("Obesity", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_stroke           = grepl("Stroke", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_chronicKidney    = grepl("Chronic kidney", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_certainCancers   = grepl("Certain cancers", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_IBD              = grepl("Inflammatory bowel", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_IBS              = grepl("Irritable bowel", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_fattyLiver       = grepl("Fatty liver", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_vascularDementia = grepl("Vascular", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_improvedConditions_depression       = grepl("Depression", Pbdiet_improvedConditions) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         diet_personallyTriedWFPB = factor(diet_personallyTriedWFPB, levels = c("No;", "Yes;"), labels = c("no", "yes")),
         diet_personallyTriedWFPB_yes_eliminatedDairy   = case_when(diet_personallyTriedWFPB == "no" ~ NA_character_,
                                                                        TRUE ~ grepl("Dairy", diet_personallyTriedWFPB_yes_eliminatedFoods) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_yes_eliminatedEggs    = case_when(diet_personallyTriedWFPB == "no" ~ NA_character_,
                                                                        TRUE ~ grepl("Eggs", diet_personallyTriedWFPB_yes_eliminatedFoods) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_yes_eliminatedFish    = case_when(diet_personallyTriedWFPB == "no" ~ NA_character_,
                                                                        TRUE ~ grepl("Fish", diet_personallyTriedWFPB_yes_eliminatedFoods) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_yes_eliminatedHoney   = case_when(diet_personallyTriedWFPB == "no" ~ NA_character_,
                                                                        TRUE ~ grepl("Honey", diet_personallyTriedWFPB_yes_eliminatedFoods) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_yes_eliminatedPoultry = case_when(diet_personallyTriedWFPB == "no" ~ NA_character_,
                                                                        TRUE ~ grepl("Poultry", diet_personallyTriedWFPB_yes_eliminatedFoods) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_yes_eliminatedRedMeat = case_when(diet_personallyTriedWFPB == "no" ~ NA_character_,
                                                                        TRUE ~ grepl("Red meat", diet_personallyTriedWFPB_yes_eliminatedFoods) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_yes_eliminatedNothing = case_when(diet_personallyTriedWFPB == "no" ~ NA_character_,
                                                                        TRUE ~ grepl("Not applicable", diet_personallyTriedWFPB_yes_eliminatedFoods) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_no_barriersPBdiet         = case_when(diet_personallyTriedWFPB == "yes" ~ NA_character_,
                                                                        TRUE ~ grepl("plant-based diet", diet_personallyTriedWFPB_no_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_no_barriersHealthy        = case_when(diet_personallyTriedWFPB == "yes" ~ NA_character_,
                                                                        TRUE ~ grepl("healthy", diet_personallyTriedWFPB_no_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_no_barriersProtein        = case_when(diet_personallyTriedWFPB == "yes" ~ NA_character_,
                                                                        TRUE ~ grepl("protein", diet_personallyTriedWFPB_no_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_no_barriersDairy          = case_when(diet_personallyTriedWFPB == "yes" ~ NA_character_,
                                                                        TRUE ~ grepl("dairy", diet_personallyTriedWFPB_no_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_no_barriersEggs           = case_when(diet_personallyTriedWFPB == "yes" ~ NA_character_,
                                                                        TRUE ~ grepl("eggs", diet_personallyTriedWFPB_no_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_no_barriersCost           = case_when(diet_personallyTriedWFPB == "yes" ~ NA_character_,
                                                                        TRUE ~ grepl("cost", diet_personallyTriedWFPB_no_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_no_barriersMicronutrients = case_when(diet_personallyTriedWFPB == "yes" ~ NA_character_,
                                                                        TRUE ~ grepl("micronutrient", diet_personallyTriedWFPB_no_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_no_barriersCultural       = case_when(diet_personallyTriedWFPB == "yes" ~ NA_character_,
                                                                        TRUE ~ grepl("cultural", diet_personallyTriedWFPB_no_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         diet_personallyTriedWFPB_no_barriersNothing        = case_when(diet_personallyTriedWFPB == "yes" ~ NA_character_,
                                                                        TRUE ~ grepl("Not applicable", diet_personallyTriedWFPB_no_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes"))),
         WFPB_isSustainableLongTerm  = factor(WFPB_isSustainableLongTerm, levels = c("Strongly disagree", "Disagree", "Not sure", "Agree", "Strongly agree")),
         WFPBrecommendation_howOften = factor(WFPBrecommendation_howOften, levels = c("Not applicable", "Never", "Rarely", "Sometimes", "Often", "Always")),
         patientsGoWFPB_barriersInterest    = grepl("interest", patientsGoWFPB_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         patientsGoWFPB_barriersDifficulty  = grepl("difficulty", patientsGoWFPB_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         patientsGoWFPB_barriersBeliefs     = grepl("beliefs", patientsGoWFPB_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         patientsGoWFPB_barriersProtein     = grepl("protein", patientsGoWFPB_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         patientsGoWFPB_barriersFoodGroups  = grepl("food groups", patientsGoWFPB_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         patientsGoWFPB_barriersCost        = grepl("cost", patientsGoWFPB_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         patientsGoWFPB_barriersKnowledge   = grepl("knowledge", patientsGoWFPB_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         patientsGoWFPB_barriersFoodOptions = grepl("Food options", patientsGoWFPB_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         patientsGoWFPB_barriersMealPrep    = grepl("Preparing meals", patientsGoWFPB_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         patientsGoWFPB_barriersEthnic      = grepl("ethnic", patientsGoWFPB_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         patientsGoWFPB_barriersNothing     = grepl("Not applicable", patientsGoWFPB_barriers) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         education_WFPBeducation     = factor(education_WFPBeducation,
                                              levels = c("Did not receive any training or resources on this topic", "Very poor", "Poor", "Good", "Very good", "Excellent"),
                                              labels = c("No training", "Very poor", "Poor", "Good", "Very good", "Excellent")),
         education_WFPBeducation_num = as.numeric(education_WFPBeducation) - 1,
         WFPBcounselling_confidence  = factor(WFPBcounselling_confidence,
                                              levels = c("Not applicable", "Not confident", "Slightly confident", "Somewhat confident",
                                                                                     "Confident as I have done it personally, not because of my training", "Fairly confident", "Completely confident"),
                                              labels = c("Not applicable", "Not confident", "Slightly confident", "Somewhat confident", "Confident", "Fairly confident", "Completely confident")),
         T2Diab_dietNational    = grepl("National", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         T2Diab_dietLowCarbs    = grepl("Low carbohydrate", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         T2Diab_dietKeto        = grepl("Ketogenic", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         T2Diab_dietMedit       = grepl("Mediterranean", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         T2Diab_dietDASH        = grepl("DASH", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         T2Diab_dietVegetarian  = grepl("Vegetarian", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         T2Diab_dietVegan       = grepl("Vegan", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         T2Diab_dietWFPB        = grepl("Whole food plant-based", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         T2Diab_dietHighProtein = grepl("High protein", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         T2Diab_dietTDR         = grepl("TDR", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         T2Diab_dietPDR         = grepl("PDR", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         T2Diab_dietLowEnergy   = grepl("Low energy", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         T2Diab_dietLowFat      = grepl("Low fat", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         T2Diab_dietNA          = grepl("Not applicable", T2Diab_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietNational    = grepl("National", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietLowCarbs    = grepl("Low carbohydrate", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietKeto        = grepl("Ketogenic", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietMedit       = grepl("Mediterranean", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietDASH        = grepl("DASH", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietHEART       = grepl("HEART", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietVegetarian  = grepl("Vegetarian", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietVegan       = grepl("Vegan", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietWFPB        = grepl("Whole food plant-based", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietHighProtein = grepl("High protein", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietTDR         = grepl("TDR", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietPDR         = grepl("PDR", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietLowEnergy   = grepl("Low energy", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietLowFat      = grepl("Low fat", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cardiovascDisease_dietNA          = grepl("Not applicable", cardiovascDisease_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietNational    = grepl("National", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietLowCarbs    = grepl("Low carbohydrate", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietKeto        = grepl("Ketogenic", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietMedit       = grepl("Mediterranean", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietDASH        = grepl("DASH", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietVegetarian  = grepl("Vegetarian", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietVegan       = grepl("Vegan", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietWFPB        = grepl("Whole food plant-based", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietHighProtein = grepl("High protein", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietTDR         = grepl("TDR", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietPDR         = grepl("PDR", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietLowEnergy   = grepl("Low energy", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietLowFat      = grepl("Low fat", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietPortfolio   = grepl("Portfolio", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         weightLoss_dietNA          = grepl("Not applicable", weightLoss_diet) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cancerPrev_guidelinesNDG      = grepl("National Dietary Guidelines", cancerPrev_dietaryGuidelines) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cancerPrev_guidelinesWCRF     = grepl("WCRF", cancerPrev_dietaryGuidelines) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cancerPrev_guidelinesWHO      = grepl("WHO", cancerPrev_dietaryGuidelines) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cancerPrev_guidelinesAmerican = grepl("American", cancerPrev_dietaryGuidelines) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cancerPrev_guidelinesESPEN    = grepl("ESPEN", cancerPrev_dietaryGuidelines) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cancerPrev_guidelinesCancerUK = grepl("Cancer Research UK", cancerPrev_dietaryGuidelines) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cancerPrev_guidelinesIrish    = grepl("Irish", cancerPrev_dietaryGuidelines) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cancerPrev_guidelinesCanadian = grepl("Canadian", cancerPrev_dietaryGuidelines) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cancerPrev_guidelinesNCI      = grepl("National Cancer Institute", cancerPrev_dietaryGuidelines) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         cancerPrev_guidelinesNA       = grepl("Not applicable", cancerPrev_dietaryGuidelines) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         WFPB_financiallyRealistic_forLowIncome  = factor(WFPB_financiallyRealistic_forLowIncome, levels = c("Strongly disagree", "Disagree", "Not sure", "Agree", "Strongly agree")),
         Pbdiet_enoughProperEducationalResources = case_when(Pbdiet_enoughProperEducationalResources == "Could be better but there are some" ~ "Agree",
                                                             Pbdiet_enoughProperEducationalResources == "Feel more  + more are becoming available e.g. RNG plant based" ~ "Agree",
                                                             Pbdiet_enoughProperEducationalResources == "More UK based than Irish"           ~ "Not sure",
                                                             TRUE ~ Pbdiet_enoughProperEducationalResources),
         Pbdiet_enoughProperEducationalResources = factor(Pbdiet_enoughProperEducationalResources, levels = c("Strongly disagree", "Disagree", "Not sure", "Agree", "Strongly agree")),
         withinHealthcare_PbdietShouldBeOption   = case_when(grepl("agree", tolower(withinHealthcare_PbdietShouldBeOption))    ~ withinHealthcare_PbdietShouldBeOption,
                                                             grepl("expensive", withinHealthcare_PbdietShouldBeOption)         ~ "Disagree",
                                                             grepl("IF the patient", withinHealthcare_PbdietShouldBeOption)    ~ "Agree",
                                                             withinHealthcare_PbdietShouldBeOption == "it should be an option" ~ "Agree",
                                                             grepl("ADDITIONAL option", withinHealthcare_PbdietShouldBeOption) ~ "Agree",
                                                             grepl("Provision should", withinHealthcare_PbdietShouldBeOption)  ~ "Agree",
                                                             TRUE ~ "Not sure"),
         withinHealthcare_PbdietShouldBeOption   = factor(withinHealthcare_PbdietShouldBeOption, levels = c("Strongly disagree", "Disagree", "Not sure", "Agree", "Strongly agree")),
         workplace_supportForWFPBadvocacy = factor(workplace_supportForWFPBadvocacy, levels = c("Strongly unsupported", "Somewhat unsupported", "Neutral", "Somewhat supported", "Strongly supported")),
         Pbdiet_increasedRiskForMalnutrition = factor(Pbdiet_increasedRiskForMalnutrition, levels = c("Strongly disagree", "Disagree", "Not sure", "Agree", "Strongly agree")),
         Pbdiet_increasedRiskEatingDisorder = case_when(grepl("agree", tolower(Pbdiet_increasedRiskEatingDisorder))       ~ Pbdiet_increasedRiskEatingDisorder,
                                                        grepl("often connected", Pbdiet_increasedRiskEatingDisorder)      ~ "Disagree",
                                                        Pbdiet_increasedRiskEatingDisorder == "Some degree"               ~ "Agree",
                                                        grepl("Not causal", Pbdiet_increasedRiskEatingDisorder)           ~ "Disagree",
                                                        grepl("causal factor", Pbdiet_increasedRiskEatingDisorder)        ~ "Disagree",
                                                        grepl("No ", Pbdiet_increasedRiskEatingDisorder)                  ~ "Disagree",
                                                        grepl("A restrictive", Pbdiet_increasedRiskEatingDisorder)        ~ "Agree",
                                                        grepl("further developing", Pbdiet_increasedRiskEatingDisorder)   ~ "Disagree",
                                                        grepl("Doesn't increase", Pbdiet_increasedRiskEatingDisorder)     ~ "Disagree",
                                                        grepl("The dietary approach", Pbdiet_increasedRiskEatingDisorder) ~ "Disagree",
                                                        grepl("The diet itself", Pbdiet_increasedRiskEatingDisorder)      ~ "Disagree",
                                                        grepl("It does not", Pbdiet_increasedRiskEatingDisorder)          ~ "Disagree",
                                                        grepl("It is probably", Pbdiet_increasedRiskEatingDisorder)       ~ "Disagree",
                                                        grepl("People with food", Pbdiet_increasedRiskEatingDisorder)     ~ "Disagree",
                                                        grepl("It can be associated", Pbdiet_increasedRiskEatingDisorder) ~ "Disagree",
                                                        grepl("may be an indicator", Pbdiet_increasedRiskEatingDisorder)  ~ "Disagree",
                                                        grepl("I don't think it", Pbdiet_increasedRiskEatingDisorder)     ~ "Disagree",
                                                        grepl("it doesn't cause", Pbdiet_increasedRiskEatingDisorder)     ~ "Disagree",
                                                        grepl("It's not a direct", Pbdiet_increasedRiskEatingDisorder)    ~ "Disagree",
                                                        grepl("Can be used as", Pbdiet_increasedRiskEatingDisorder)       ~ "Disagree",
                                                        grepl("I know there is an", Pbdiet_increasedRiskEatingDisorder)   ~ "Disagree",
                                                        TRUE ~ "Not sure"),
         Pbdiet_increasedRiskEatingDisorder = factor(Pbdiet_increasedRiskEatingDisorder, levels = c("Strongly disagree", "Disagree", "Not sure", "Agree", "Strongly agree")),
         Pbdiet_primaryMotivationClientsWeightLoss            = grepl("Weight loss", Pbdiet_primaryMotivationClients) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_primaryMotivationClientsCardiometabolicHealth = grepl("Cardiometabolic health", Pbdiet_primaryMotivationClients) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_primaryMotivationClientsOverallHealth         = grepl("overall health", Pbdiet_primaryMotivationClients) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_primaryMotivationClientsGlycaemicIndex        = grepl("glycaemic", Pbdiet_primaryMotivationClients) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_primaryMotivationClientsEnvironmental         = grepl("Environmental", Pbdiet_primaryMotivationClients) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_primaryMotivationClientsEthics                = grepl("Ethics", Pbdiet_primaryMotivationClients) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         Pbdiet_primaryMotivationClientsNothing               = grepl("Not applicable", Pbdiet_primaryMotivationClients) %>% factor(levels = c(FALSE, TRUE), labels = c("no", "yes")),
         ChronicConditions_longTermWFPBadherence   = case_when(grepl("agree", tolower(ChronicConditions_longTermWFPBadherence))      ~ withinHealthcare_PbdietShouldBeOption,
                                                              grepl("support", ChronicConditions_longTermWFPBadherence)             ~ "Agree",
                                                              grepl("If their motivation", ChronicConditions_longTermWFPBadherence) ~ "Agree",
                                                              grepl("Only if it was", ChronicConditions_longTermWFPBadherence)      ~ "Agree",
                                                              grepl("Not likely if", ChronicConditions_longTermWFPBadherence)       ~ "Agree",
                                                              grepl("Not appropriate", ChronicConditions_longTermWFPBadherence)     ~ "Disagree",
                                                              TRUE ~ "Not sure"),
         ChronicConditions_longTermWFPBadherence = factor(ChronicConditions_longTermWFPBadherence, levels = c("Strongly disagree", "Disagree", "Not sure", "Agree", "Strongly agree"))
         )
Code
# define a vector of the subspecialty areas that we focus on
subspecialties_varVector <- c("weight management" = "dieteticsSubspecialty_obesityOrWeightManagement",
                              "gastroenterology"  = "dieteticsSubspecialty_gastroenterology",
                              "diabetes"          = "dieteticsSubspecialty_diabetes",
                              "paediatrics"       = "dieteticsSubspecialty_paediatrics",
                              "elderly care"      = "dieteticsSubspecialty_elderly",
                              "oncology"          = "dieteticsSubspecialty_oncology",
                              "eating disorders"  = "dieteticsSubspecialty_psychEating")
Code
# divide age and years of practice by 10, to facilitate the interpretation of regression coefficients
dat <- dat %>% 
  mutate(age_10yearUnits                  = age / 10,
         dietitian_forHowLong_10yearUnits = dietitian_forHowLong / 10)

Description baseline demographics

For the diet variable, ‘Flexitarian’ and ‘Pescetarian’ are combined into one group.

Code
dat_descr <- dat %>% 
  mutate(diet = case_when(diet_personal %in% c("Flexitarian", "Pescetarian") ~ "Semi-vegetarian",
                          TRUE ~ diet_personal),
         diet = factor(diet,
                       levels = c("Omnivore", "Mediterranean", "Semi-vegetarian",
                                  "Vegetarian", "WFPB", "Vegan")))
Code
# helper function to nicely print a table
print_table <- function(tab) {
  tab %>% 
    as.data.frame() %>% 
    rename(frequency = Freq) %>% 
    mutate(rel_frequency = frequency / sum(frequency),
           rel_frequency = paste0(100 * round(rel_frequency, 2), "%")) %>% 
    kable() %>% 
    kable_styling()
}
Code
# helper function to paste the matrices of absolute and relative frequencies together
paste_tabs <- function(tab_abs, tab_rel) {
  
  rel_vec <- paste0("(", 100 * round(tab_rel, 2), "%)")
  rel_vec <- gsub(rel_vec, pattern = "NaN%", replacement = "-")
  
  paste(tab_abs, rel_vec) %>%
    matrix(nrow     = nrow(tab_abs),
           ncol     = ncol(tab_abs),
           byrow    = FALSE,
           dimnames = list(row.names(tab_abs),
                           colnames(tab_abs)))
}

Diet

Overall diet distribution

Code
# diet frequency table
table(dat_descr$diet, useNA = "always") %>% print_table()
Var1 frequency rel_frequency
Omnivore 103 31%
Mediterranean 56 17%
Semi-vegetarian 102 30%
Vegetarian 30 9%
WFPB 20 6%
Vegan 24 7%
NA 0 0%

Gender

Gender overall

The one NA value is a ‘prefer not to say’ answer.

Code
table(dat_descr$gender, useNA = "always") %>% print_table()
Var1 frequency rel_frequency
Male 10 3%
Female 324 97%
NA 1 0%

Gender by diet

Code
tab_abs <- table(dat_descr$gender, dat_descr$diet, useNA = "always")
tab_rel <- tab_abs %>% prop.table(margin = 2)

paste_tabs(tab_abs, tab_rel) %>% kable() %>% kable_styling()
Omnivore Mediterranean Semi-vegetarian Vegetarian WFPB Vegan NA
Male 5 (5%) 3 (5%) 1 (1%) 0 (0%) 1 (5%) 0 (0%) 0 (-)
Female 97 (94%) 53 (95%) 101 (99%) 30 (100%) 19 (95%) 24 (100%) 0 (-)
NA 1 (1%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (-)

Age

Age overall

The one NA value is an implausible ‘age 2 years’ answer.

Code
dat_descr <- dat_descr %>% 
  mutate(age_newCat = case_when(is.na(age) ~ NA_character_,
                                age < 30   ~ "20-29",
                                age < 40   ~ "30-39",
                                age < 50   ~ "40-49",
                                TRUE       ~ ">= 50"),
         age_newCat = factor(age_newCat,
                             levels = c("20-29","30-39","40-49",">= 50")))
Code
table(dat_descr$age_newCat, useNA = "always") %>% print_table()
Var1 frequency rel_frequency
20-29 73 22%
30-39 108 32%
40-49 87 26%
>= 50 66 20%
NA 1 0%
Code
median_age <- median(dat_descr$age, na.rm = TRUE)

Median age overall: 38

Age by diet

Code
tab_abs <- table(dat_descr$age_newCat, dat_descr$diet, useNA = "always")
tab_rel <- tab_abs %>% prop.table(margin = 2)

paste_tabs(tab_abs, tab_rel) %>% kable() %>% kable_styling()
Omnivore Mediterranean Semi-vegetarian Vegetarian WFPB Vegan NA
20-29 20 (19%) 9 (16%) 24 (24%) 9 (30%) 4 (20%) 7 (29%) 0 (-)
30-39 38 (37%) 15 (27%) 33 (32%) 7 (23%) 8 (40%) 7 (29%) 0 (-)
40-49 23 (22%) 14 (25%) 28 (27%) 8 (27%) 4 (20%) 10 (42%) 0 (-)
>= 50 21 (20%) 18 (32%) 17 (17%) 6 (20%) 4 (20%) 0 (0%) 0 (-)
NA 1 (1%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (-)
Code
dat_descr %>% 
  group_by(diet) %>% 
  summarize(median_age = median(age, na.rm = TRUE)) %>% 
  kable() %>% 
  kable_styling()
diet median_age
Omnivore 38.0
Mediterranean 42.5
Semi-vegetarian 36.5
Vegetarian 36.0
WFPB 37.0
Vegan 35.0

Education

Education overall

Code
table(dat_descr$education, useNA = "always") %>% print_table()
Var1 frequency rel_frequency
Undergraduate degree 142 42%
Postgraduate degree 172 51%
PhD 20 6%
NA 1 0%

Education by diet

Code
tab_abs <- table(dat_descr$education, dat_descr$diet, useNA = "always")
tab_rel <- tab_abs %>% prop.table(margin = 2)

paste_tabs(tab_abs, tab_rel) %>% kable() %>% kable_styling()
Omnivore Mediterranean Semi-vegetarian Vegetarian WFPB Vegan NA
Undergraduate degree 45 (44%) 22 (39%) 45 (44%) 11 (37%) 9 (45%) 10 (42%) 0 (-)
Postgraduate degree 48 (47%) 34 (61%) 49 (48%) 17 (57%) 11 (55%) 13 (54%) 0 (-)
PhD 10 (10%) 0 (0%) 8 (8%) 1 (3%) 0 (0%) 1 (4%) 0 (-)
NA 0 (0%) 0 (0%) 0 (0%) 1 (3%) 0 (0%) 0 (0%) 0 (-)

Years of practice

Years of practice overall

Code
table(dat_descr$dietitian_forHowLong_cat, useNA = "always") %>% print_table()
Var1 frequency rel_frequency
up to 3 years 75 22%
4 to 6 years 63 19%
7 to 9 years 34 10%
10 to 14 years 52 16%
15 to 19 years 34 10%
20 to 24 years 29 9%
25+ years 48 14%
NA 0 0%
Code
median_yop <- median(dat_descr$dietitian_forHowLong, na.rm = TRUE)

Median years of practice overall: 9

Years of practice by diet

Code
tab_abs <- table(dat_descr$dietitian_forHowLong_cat, dat_descr$diet, useNA = "always")
tab_rel <- tab_abs %>% prop.table(margin = 2)

paste_tabs(tab_abs, tab_rel) %>% kable() %>% kable_styling()
Omnivore Mediterranean Semi-vegetarian Vegetarian WFPB Vegan NA
up to 3 years 25 (24%) 11 (20%) 21 (21%) 7 (23%) 2 (10%) 9 (38%) 0 (-)
4 to 6 years 18 (17%) 8 (14%) 19 (19%) 6 (20%) 8 (40%) 4 (17%) 0 (-)
7 to 9 years 12 (12%) 5 (9%) 7 (7%) 3 (10%) 2 (10%) 5 (21%) 0 (-)
10 to 14 years 16 (16%) 10 (18%) 18 (18%) 0 (0%) 4 (20%) 4 (17%) 0 (-)
15 to 19 years 8 (8%) 6 (11%) 13 (13%) 4 (13%) 1 (5%) 2 (8%) 0 (-)
20 to 24 years 9 (9%) 7 (12%) 10 (10%) 3 (10%) 0 (0%) 0 (0%) 0 (-)
25+ years 15 (15%) 9 (16%) 14 (14%) 7 (23%) 3 (15%) 0 (0%) 0 (-)
NA 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (-)
Code
dat_descr %>% 
  group_by(diet) %>% 
  summarize(median_years = median(dietitian_forHowLong, na.rm = TRUE)) %>% 
  kable() %>% 
  kable_styling()
diet median_years
Omnivore 8.0
Mediterranean 11.5
Semi-vegetarian 10.5
Vegetarian 8.0
WFPB 6.5
Vegan 5.5

Area of work

Area of work overall

Code
tab_workArea <- dat_descr %>% 
  summarize(n_hospital            = sum(dieteticsArea_hospital == "yes"),
            n_primCareORcommunity = sum((dieteticsArea_primaryCare == "yes") |
                                          (dieteticsArea_community == "yes")),
            n_private             = sum(dieteticsArea_privatePractice == "yes"),
            n_academia            = sum(dieteticsArea_academia == "yes"),
            n_publicHealth        = sum(dieteticsArea_publicHealth == "yes"))

tab_workArea %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(frequency = ".") %>% 
  mutate(rel_frequency = frequency / sum(frequency),
         rel_frequency = paste0(100 * round(rel_frequency, 2), "%")) %>% 
  kable() %>% 
  kable_styling()
frequency rel_frequency
n_hospital 148 40%
n_primCareORcommunity 135 36%
n_private 42 11%
n_academia 25 7%
n_publicHealth 22 6%

Area of work by diet

Code
tab_workArea <- dat_descr %>% 
  group_by(diet) %>% 
  summarize(n_hospital            = sum(dieteticsArea_hospital == "yes"),
            n_primCareORcommunity = sum((dieteticsArea_primaryCare == "yes") |
                                          (dieteticsArea_community == "yes")),
            n_private             = sum(dieteticsArea_privatePractice == "yes"),
            n_academia            = sum(dieteticsArea_academia == "yes"),
            n_publicHealth        = sum(dieteticsArea_publicHealth == "yes"))
tab_abs <- unlist(tab_workArea %>% select(-diet), use.names = FALSE) %>% 
  matrix(nrow     = ncol(tab_workArea) - 1,
         ncol     = nrow(tab_workArea),
         byrow    = TRUE,
         dimnames = list(colnames(tab_workArea)[-1],
                         tab_workArea$diet))
tab_rel <- tab_abs %>% prop.table(margin = 2)

paste_tabs(tab_abs, tab_rel) %>% kable() %>% kable_styling()
Omnivore Mediterranean Semi-vegetarian Vegetarian WFPB Vegan
n_hospital 47 (40%) 18 (32%) 54 (48%) 9 (27%) 8 (35%) 12 (41%)
n_primCareORcommunity 44 (37%) 29 (52%) 35 (31%) 10 (30%) 8 (35%) 9 (31%)
n_private 7 (6%) 5 (9%) 11 (10%) 7 (21%) 6 (26%) 6 (21%)
n_academia 12 (10%) 2 (4%) 7 (6%) 3 (9%) 0 (0%) 1 (3%)
n_publicHealth 8 (7%) 2 (4%) 6 (5%) 4 (12%) 1 (4%) 1 (3%)

Area of specialty

Area of specialty overall

Code
tab_subsp <- dat_descr %>% 
  summarize(n_obesity         = sum(dieteticsSubspecialty_obesityOrWeightManagement == "yes"),
            n_diabetes        = sum(dieteticsSubspecialty_diabetes == "yes"),
            n_gastroen        = sum(dieteticsSubspecialty_gastroenterology == "yes"),
            n_paediatr        = sum(dieteticsSubspecialty_paediatrics == "yes"),
            n_elderly         = sum(dieteticsSubspecialty_elderly == "yes"),
            n_oncology        = sum(dieteticsSubspecialty_oncology == "yes"),
            n_eatingDisorders = sum(dieteticsSubspecialty_psychEating == "yes"))

tab_subsp %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(frequency = ".") %>% 
  mutate(rel_frequency = frequency / sum(frequency),
         rel_frequency = paste0(100 * round(rel_frequency, 2), "%")) %>% 
  kable() %>% 
  kable_styling()
frequency rel_frequency
n_obesity 83 22%
n_diabetes 59 16%
n_gastroen 57 15%
n_paediatr 54 14%
n_elderly 51 14%
n_oncology 36 10%
n_eatingDisorders 35 9%

Area of specialty by diet

Code
tab_subsp <- dat_descr %>% 
  group_by(diet) %>% 
  summarize(n_obesity         = sum(dieteticsSubspecialty_obesityOrWeightManagement == "yes"),
            n_diabetes        = sum(dieteticsSubspecialty_diabetes == "yes"),
            n_gastroen        = sum(dieteticsSubspecialty_gastroenterology == "yes"),
            n_paediatr        = sum(dieteticsSubspecialty_paediatrics == "yes"),
            n_elderly         = sum(dieteticsSubspecialty_elderly == "yes"),
            n_oncology        = sum(dieteticsSubspecialty_oncology == "yes"),
            n_eatingDisorders = sum(dieteticsSubspecialty_psychEating == "yes"))
tab_abs <- unlist(tab_subsp %>% select(-diet), use.names = FALSE) %>% 
  matrix(nrow     = ncol(tab_subsp) - 1,
         ncol     = nrow(tab_subsp),
         byrow    = TRUE,
         dimnames = list(colnames(tab_subsp)[-1],
                         tab_subsp$diet))
tab_rel <- tab_abs %>% prop.table(margin = 2)

paste_tabs(tab_abs, tab_rel) %>% kable() %>% kable_styling()
Omnivore Mediterranean Semi-vegetarian Vegetarian WFPB Vegan
n_obesity 20 (17%) 18 (26%) 22 (22%) 8 (29%) 8 (28%) 7 (24%)
n_diabetes 16 (14%) 8 (12%) 19 (19%) 6 (21%) 8 (28%) 2 (7%)
n_gastroen 20 (17%) 7 (10%) 16 (16%) 6 (21%) 4 (14%) 4 (14%)
n_paediatr 18 (15%) 8 (12%) 17 (17%) 5 (18%) 1 (3%) 5 (17%)
n_elderly 18 (15%) 9 (13%) 11 (11%) 2 (7%) 5 (17%) 6 (21%)
n_oncology 15 (13%) 10 (14%) 9 (9%) 0 (0%) 1 (3%) 1 (3%)
n_eatingDisorders 11 (9%) 9 (13%) 8 (8%) 1 (4%) 2 (7%) 4 (14%)

Tests baseline demographics

For all of the following baseline demographics a hypothesis test is calculated to test for potential differences between dietitians’ diets, and the respective p-value is reported. For metric variables (age and years of practice) an ANOVA is calculated, for categorical variables a chi-squared test of independence is calculated.
For the diet variable, ‘Flexitarian’ and ‘Pescetarian’ are combined into one group.

Code
# test calculation
m_age <- lm(age ~ diet, data = dat_descr)
p_age <- anova(m_age)$`Pr(>F)`[1]
p_edu <- chisq.test(x = dat_descr$diet,
                    y = dat_descr$education)$p.value
m_ypr <- lm(dietitian_forHowLong ~ diet, data = dat_descr)
p_ypr <- anova(m_ypr)$`Pr(>F)`[1]
# Note: both the 'area of work' and 'area of subspecialty' are multiple choice items, so we first calculate a contingency table on which we then run the test.
# Also, only selected areas are evaluated for both
tab_workArea <- dat_descr %>% 
  group_by(diet) %>% 
  summarize(n_hospital            = sum(dieteticsArea_hospital == "yes"),
            n_primCareORcommunity = sum((dieteticsArea_primaryCare == "yes") |
                                          (dieteticsArea_community == "yes")),
            n_private             = sum(dieteticsArea_privatePractice == "yes"),
            n_academia            = sum(dieteticsArea_academia == "yes"),
            n_publicHealth        = sum(dieteticsArea_publicHealth == "yes")) %>% 
  select(-diet)
p_workArea <- chisq.test(tab_workArea)$p.value
tab_subsp <- dat_descr %>% 
  group_by(diet) %>% 
  summarize(n_obesity         = sum(dieteticsSubspecialty_obesityOrWeightManagement == "yes"),
            n_diabetes        = sum(dieteticsSubspecialty_diabetes == "yes"),
            n_gastroen        = sum(dieteticsSubspecialty_gastroenterology == "yes"),
            n_paediatr        = sum(dieteticsSubspecialty_paediatrics == "yes"),
            n_elderly         = sum(dieteticsSubspecialty_elderly == "yes"),
            n_oncology        = sum(dieteticsSubspecialty_oncology == "yes"),
            n_eatingDisorders = sum(dieteticsSubspecialty_psychEating == "yes")) %>% 
  select(-diet)
p_subsp <- chisq.test(tab_subsp)$p.value
  • test between diet and sex: not tested due to too small n
  • test between diet and age: 0.0442
  • test between diet and level of education: 0.4079
  • test between diet and years of working as a dietitian: 0.088
  • test between diet and area of work: 0.056
  • test between diet and area of subspecialty: 0.5155

Regression model preparation

Exclude purely academic dietitians

6 dietitians have a purely academic area of work. These dietitians will be excluded from all following analyses since focus is on practitioning dietitians.

Code
dat <- dat %>% filter(dietetics_area != "Academia;")

Sensitivity analysis: Exclusion of new dietitians

As a sensitivity analysis we report the result figures in the supplementary materials of the publication for the case when excluding all newly practitioning dietitians (i.e., dietitians with less than 3 years of practical experience) from the estimation. In this way we can check if the inclusion or exclusion of this group influences the study results to a relevant extent.

By default, all of the following analyses are calculated without excluding this group of dietitians. If you want to exclude them and thus run the sensitivity analysis, simply uncomment the following line of code and re-run all of the following analyses.

Code
# dat <- dat %>% filter(dietitian_forHowLong >= 3)

Check for multicollinearities between subspecialties

We focus on six subspecialties in all analyses. As can be seen in the following matrix, there is no problematic multicollinearity present between them:

Code
cor_tab <- dat %>% 
  select(all_of(unname(subspecialties_varVector))) %>% 
  mutate_all(as.numeric) %>% 
  cor() %>% 
  round(2)
row.names(cor_tab) <- gsub(row.names(cor_tab), pattern = "dieteticsSubspecialty_", replacement = "")
colnames(cor_tab)  <- gsub(colnames(cor_tab),  pattern = "dieteticsSubspecialty_", replacement = "")
cor_tab %>% 
  kable() %>% 
  kable_styling()
obesityOrWeightManagement gastroenterology diabetes paediatrics elderly oncology psychEating
obesityOrWeightManagement 1.00 0.01 0.36 -0.07 -0.01 -0.06 -0.02
gastroenterology 0.01 1.00 -0.09 -0.01 0.12 0.02 0.00
diabetes 0.36 -0.09 1.00 -0.05 0.00 -0.06 -0.11
paediatrics -0.07 -0.01 -0.05 1.00 -0.10 -0.05 0.01
elderly -0.01 0.12 0.00 -0.10 1.00 -0.01 -0.04
oncology -0.06 0.02 -0.06 -0.05 -0.01 1.00 -0.06
psychEating -0.02 0.00 -0.11 0.01 -0.04 -0.06 1.00

Define helper functions

Code
# helper function for model estimation
# - dat: dataset
# - y:   character name of the response variable
# The function returns a list of the model results and the AUC value
estimate_logitModel <- function(dat, y) {
  
  # rename response variable, for easier handling
  colnames(dat)[colnames(dat) == y] <- "y"
  
  # model estimation
  model <- gam(y ~ 
                 dieteticsSubspecialty_obesityOrWeightManagement +
                 dieteticsSubspecialty_gastroenterology +
                 dieteticsSubspecialty_diabetes +
                 dieteticsSubspecialty_paediatrics +
                 dieteticsSubspecialty_elderly +
                 dieteticsSubspecialty_oncology +
                 dieteticsSubspecialty_psychEating +
                 dietitian_forHowLong_10yearUnits +
                 education +
                 education_WFPBeducation_num,
               family = binomial(link = "logit"),
               data = dat)
  
  # prepare the results table
  results_tab <- summary(model)$p.table %>% 
    as.data.frame() %>% 
    mutate(parameter = row.names(.)) %>% 
    rename(estimate = "Estimate",
           se       = "Std. Error",
           pvalue   = "Pr(>|z|)") %>% 
    mutate(CI_lower     = estimate - qnorm(.975) * se,
           CI_upper     = estimate + qnorm(.975) * se,
           exp_estimate = exp(estimate),
           exp_CIlower  = exp(CI_lower),
           exp_CIupper  = exp(CI_upper),
           pvalue       = case_when(pvalue < .0001 ~ "<.0001",
                                    TRUE           ~ as.character(round(pvalue, 4))),
           param_group  = case_when(grepl("Subspecialty", parameter)  ~ "area of specialty",
                                    grepl("forHowLong", parameter)    ~ "years of practice",
                                    grepl("WFPBeducation", parameter) ~ "WFPB education",
                                    grepl("education", parameter)     ~ "education"),
           param_group  = factor(param_group, levels = c("area of specialty", "years of practice", "education", "WFPB education")),
           parameter    = gsub(parameter, pattern = "dieteticsSubspecialty_", replacement = "")) %>% 
    select(param_group, parameter, exp_estimate, exp_CIlower, exp_CIupper, pvalue)
  row.names(results_tab) <- NULL
  
  # The model's goodness-of-fit is determined by calculating the AUC (Area Under the Curve).
  # Therefore, the dataset is first randomly split into a training dataset (80% of the data)
  # and a test dataset (20% of the data).
  # The regression model is then re-estimated on the training data,
  # and the AUC is calculated on the test data.
  
  # basis: split data into a training set (80% of the data) and a test set (20%)
  set.seed(2024)
  train_obs <- sample(1:nrow(dat), size = round(0.8*nrow(dat)), replace = FALSE)
  dat_train <- dat %>% slice(train_obs)
  dat_test  <- dat %>% slice(-train_obs)
  
  # 1) estimate the model on the training set
  model_train <- gam(y ~ 
                       dieteticsSubspecialty_obesityOrWeightManagement +
                       dieteticsSubspecialty_gastroenterology +
                       dieteticsSubspecialty_diabetes +
                       dieteticsSubspecialty_paediatrics +
                       dieteticsSubspecialty_elderly +
                       dieteticsSubspecialty_oncology +
                       dieteticsSubspecialty_psychEating +
                       dietitian_forHowLong_10yearUnits +
                       education,
                     family = binomial(link = "logit"),
                     data = dat_train)
  # 2) calculate ROC curve and the corresponding AUC value
  pred <- data.frame("outcome" = as.integer(dat_test$y) - 1,
                     "pred"    = predict(model_train, newdata = dat_test))
  gg <- ggplot(pred, aes(m = pred, d = outcome)) + geom_roc() + ggtitle("ROC curve")
  auc <- calc_auc(gg)$AUC
  
  
  # return results
  list("results_tab" = results_tab,
       "AUC"         = auc)
}
Code
# helper function to create an effect plot, based on a regression results table
plot_logitRegressionResults <- function(results_tab) {
  # effect plot
  results_tab %>% 
    filter(parameter != "(Intercept)") %>% 
    ggplot(aes(x = parameter)) +
    geom_hline(yintercept = 1, lty = 2, col = "gray50") +
    geom_pointrange(aes(y = exp_estimate, ymin = exp_CIlower, ymax = exp_CIupper, col = param_group)) +
    scale_y_continuous("Odds Ratio, on log2-transformed scale", trans = "log2") +
    facet_wrap(~ param_group, scales = "free_x", nrow = 1) +
    theme(axis.text.x     = element_text(angle = 45, hjust = 1),
          axis.title.x    = element_blank(),
          legend.position = "none")
}

Knowledge research questions

Plant-based: Diet definition

What is their understanding on what constitutes/defines a plant-based dietary pattern?

Code
dat %>% 
  select(starts_with("diet_perceptionPBdiet_")) %>% 
  pivot_longer(cols = everything()) %>% 
  group_by(name) %>% 
  mutate(sum_yes = sum(value == "yes")) %>% 
  ungroup() %>% 
  arrange(desc(sum_yes)) %>% 
  mutate(name = gsub(pattern = "diet_perceptionPBdiet_", replacement = "", name),
         name = factor(name, levels = rev(unique(name)))) %>% 
  ggplot(aes(y = name, fill = value)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("lightgray", unname(col_vector["knowledge"]))) +
  scale_x_continuous(labels = scales::label_percent()) +
  xlab("Agreement [%]") +
  ggtitle("Perception of what is a plant-based diet") +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

Plant-based: Risk-reduction regarding specific diseases

What is their understanding and knowledge on the role of plant-based diets in the context of reducing the risk of different non-communicable-related diseases?

Code
dat %>% 
  select(starts_with("Pbdiet_improvedConditions_")) %>% 
  pivot_longer(cols = everything()) %>% 
  group_by(name) %>% 
  mutate(sum_yes = sum(value == "yes")) %>% 
  ungroup() %>% 
  arrange(desc(sum_yes)) %>% 
  mutate(name = gsub(pattern = "Pbdiet_improvedConditions_", replacement = "", name),
         name = factor(name, levels = rev(unique(name)))) %>% 
  ggplot(aes(y = name, fill = value)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("lightgray", unname(col_vector["knowledge"]))) +
  scale_x_continuous(labels = scales::label_percent()) +
  xlab("Agreement [%]") +
  ggtitle("Improved conditions through plant-based diet") +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

Code
# T2Diab_diet
plot_dat_T2Diab <- dat %>% 
  select(starts_with("T2Diab_diet")) %>% 
  select(-T2Diab_diet) %>% 
  pivot_longer(cols = everything()) %>% 
  group_by(name) %>% 
  mutate(sum_yes = sum(value == "yes")) %>% 
  ungroup() %>% 
  arrange(desc(sum_yes)) %>% 
  mutate(name = gsub(pattern = "T2Diab_diet", replacement = "", name),
         name = factor(name, levels = rev(unique(name))))
gg1 <- plot_dat_T2Diab %>% 
  ggplot(aes(y = name, fill = value)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("lightgray", unname(col_vector["knowledge"]))) +
  scale_x_continuous(labels = scales::label_percent()) +
  xlab("Agreement [%]") +
  ggtitle("Recommended diet\nfor type-2 diabetes") +
  theme(axis.title.y = element_blank(),
        legend.title = element_blank())

# cardiovascDisease_diet
plot_dat_cardiovascDisease <- dat %>% 
  select(starts_with("cardiovascDisease_diet")) %>% 
  select(-cardiovascDisease_diet) %>% 
  pivot_longer(cols = everything()) %>% 
  group_by(name) %>% 
  mutate(sum_yes = sum(value == "yes")) %>% 
  ungroup() %>% 
  arrange(desc(sum_yes)) %>% 
  mutate(name = gsub(pattern = "cardiovascDisease_diet", replacement = "", name),
         name = factor(name, levels = rev(unique(name))))
gg2 <- plot_dat_cardiovascDisease %>% 
  ggplot(aes(y = name, fill = value)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("lightgray", unname(col_vector["knowledge"]))) +
  scale_x_continuous(labels = scales::label_percent()) +
  xlab("Agreement [%]") +
  ggtitle("Recommended diet\nfor cardiovascular disease") +
  theme(axis.title.y = element_blank(),
        legend.title = element_blank())

# weightLoss_diet
plot_dat_weightLoss <- dat %>% 
  select(starts_with("weightLoss_diet")) %>% 
  select(-weightLoss_diet) %>% 
  pivot_longer(cols = everything()) %>% 
  group_by(name) %>% 
  mutate(sum_yes = sum(value == "yes")) %>% 
  ungroup() %>% 
  arrange(desc(sum_yes)) %>% 
  mutate(name = gsub(pattern = "weightLoss_diet", replacement = "", name),
         name = factor(name, levels = rev(unique(name))))
gg3 <- plot_dat_weightLoss %>% 
  ggplot(aes(y = name, fill = value)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("lightgray", unname(col_vector["knowledge"]))) +
  scale_x_continuous(labels = scales::label_percent()) +
  xlab("Agreement [%]") +
  ggtitle("Recommended diet\nfor weight loss") +
  theme(axis.title.y = element_blank(),
        legend.title = element_blank())

# joint graphic
(gg1 + gg2 + gg3) + patchwork::plot_layout(guides = "collect")

Univarite description

by subspecialty, age, years of practice.

Code
# prepare the datasets
plot_dat_list <- lapply(1:length(subspecialties_varVector), function(i) {
  
  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # filter the dataset for dietitians with the specific subspecialty
  dat_subsp <- dat
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  # T2Diab_diet
  dat_T2Diab <- dat_subsp %>% 
    select(starts_with("T2Diab_diet")) %>% 
    select(-T2Diab_diet) %>% 
    pivot_longer(cols = everything(), names_to = "diet") %>% 
    group_by(diet) %>% 
    summarize(share_yes = sum(value == "yes") / n()) %>% 
    arrange(desc(share_yes)) %>% 
    mutate(diet         = gsub(pattern = "T2Diab_diet", replacement = "", x = diet),
           diet         = factor(diet, levels = diet),
           subspecialty = subsp_label,
           share_label  = paste0(round(100 * share_yes), "%")) %>% 
    filter(diet != "NA")

  # cardiovascDisease_diet
  dat_cardiovascDisease <- dat_subsp %>% 
    select(starts_with("cardiovascDisease_diet")) %>% 
    select(-cardiovascDisease_diet) %>% 
    pivot_longer(cols = everything(), names_to = "diet") %>% 
    group_by(diet) %>% 
    summarize(share_yes = sum(value == "yes") / n()) %>% 
    arrange(desc(share_yes)) %>% 
    mutate(diet         = gsub(pattern = "cardiovascDisease_diet", replacement = "", x = diet),
           diet         = factor(diet, levels = diet),
           subspecialty = subsp_label,
           share_label  = paste0(round(100 * share_yes), "%")) %>%  
    filter(diet != "NA")
  
  # weightLoss_diet
  dat_weightLoss <- dat_subsp %>% 
    select(starts_with("weightLoss_diet")) %>% 
    select(-weightLoss_diet) %>% 
    pivot_longer(cols = everything(), names_to = "diet") %>% 
    group_by(diet) %>% 
    summarize(share_yes = sum(value == "yes") / n()) %>% 
    arrange(desc(share_yes)) %>% 
    mutate(diet         = gsub(pattern = "weightLoss_diet", replacement = "", x = diet),
           diet         = factor(diet, levels = diet),
           subspecialty = subsp_label,
           share_label  = paste0(round(100 * share_yes), "%")) %>%  
    filter(diet != "NA")
  
  list(dat_T2Diab            = dat_T2Diab,
       dat_cardiovascDisease = dat_cardiovascDisease,
       dat_weightLoss        = dat_weightLoss) %>% 
    return()
})

# merge prepared datasets
plot_dat_T2Diab_sub <- lapply(plot_dat_list, function(x) { x$dat_T2Diab }) %>%
  dplyr::bind_rows() %>% 
  mutate(diet         = factor(diet,         levels = levels(plot_dat_T2Diab$name)),
         subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)))
plot_dat_cardiovascDisease_sub <- lapply(plot_dat_list, function(x) { x$dat_cardiovascDisease }) %>%
  dplyr::bind_rows() %>% 
  mutate(diet         = factor(diet,         levels = levels(plot_dat_cardiovascDisease$name)),
         subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)))
plot_dat_weightLoss_sub <- lapply(plot_dat_list, function(x) { x$dat_weightLoss }) %>%
  dplyr::bind_rows() %>% 
  mutate(diet         = factor(diet,         levels = levels(plot_dat_weightLoss$name)),
         subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)))

# create heatmaps
axis_limits <- c(plot_dat_T2Diab_sub$share_yes, plot_dat_cardiovascDisease_sub$share_yes,
                 plot_dat_weightLoss_sub$share_yes) %>% 
  range()

gg1 <- plot_dat_T2Diab_sub %>% 
  ggplot(aes(x = subspecialty, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  ggtitle("Recommended diet\nfor type-2 diabetes") +
  theme(axis.text.x  = element_text(angle = 45, hjust = 1),
        axis.title.x = element_blank())

gg2 <- plot_dat_cardiovascDisease_sub %>% 
  ggplot(aes(x = subspecialty, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  ggtitle("Recommended diet\nfor cardiovascular disease") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

gg3 <- plot_dat_weightLoss_sub %>% 
  ggplot(aes(x = subspecialty, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  ggtitle("Recommended diet\nfor weight loss") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

# joint graphic
(gg1 + gg2 + gg3) +
  patchwork::plot_layout(guides = "collect") +
  patchwork::plot_annotation(title = "Recommended diets by subspecialty")

Code
# prepare the datasets
plot_dat_list <- lapply(levels(dat$age_cat), function(x) {
  
  # filter the dataset for dietitians with the specific age category
  dat_age <- dat %>% filter(age_cat == x)
  
  # T2Diab_diet
  dat_T2Diab <- dat_age %>% 
    select(starts_with("T2Diab_diet")) %>% 
    select(-T2Diab_diet) %>% 
    pivot_longer(cols = everything(), names_to = "diet") %>% 
    group_by(diet) %>% 
    summarize(share_yes = sum(value == "yes") / n()) %>% 
    arrange(desc(share_yes)) %>% 
    mutate(diet        = gsub(pattern = "T2Diab_diet", replacement = "", x = diet),
           diet        = factor(diet, levels = diet),
           age_cat     = x,
           share_label = paste0(round(100 * share_yes), "%")) %>% 
    filter(diet != "NA")

  # cardiovascDisease_diet
  dat_cardiovascDisease <- dat_age %>% 
    select(starts_with("cardiovascDisease_diet")) %>% 
    select(-cardiovascDisease_diet) %>% 
    pivot_longer(cols = everything(), names_to = "diet") %>% 
    group_by(diet) %>% 
    summarize(share_yes = sum(value == "yes") / n()) %>% 
    arrange(desc(share_yes)) %>% 
    mutate(diet        = gsub(pattern = "cardiovascDisease_diet", replacement = "", x = diet),
           diet        = factor(diet, levels = diet),
           age_cat     = x,
           share_label = paste0(round(100 * share_yes), "%")) %>%  
    filter(diet != "NA")
  
  # weightLoss_diet
  dat_weightLoss <- dat_age %>% 
    select(starts_with("weightLoss_diet")) %>% 
    select(-weightLoss_diet) %>% 
    pivot_longer(cols = everything(), names_to = "diet") %>% 
    group_by(diet) %>% 
    summarize(share_yes = sum(value == "yes") / n()) %>% 
    arrange(desc(share_yes)) %>% 
    mutate(diet        = gsub(pattern = "weightLoss_diet", replacement = "", x = diet),
           diet        = factor(diet, levels = diet),
           age_cat     = x,
           share_label = paste0(round(100 * share_yes), "%")) %>%  
    filter(diet != "NA")
  
  list(dat_T2Diab            = dat_T2Diab,
       dat_cardiovascDisease = dat_cardiovascDisease,
       dat_weightLoss        = dat_weightLoss) %>% 
    return()
})

# merge prepared datasets
plot_dat_T2Diab_age <- lapply(plot_dat_list, function(x) { x$dat_T2Diab }) %>%
  dplyr::bind_rows() %>% 
  mutate(diet    = factor(diet,    levels = levels(plot_dat_T2Diab$name)),
         age_cat = factor(age_cat, levels = levels(dat$age_cat)))
plot_dat_cardiovascDisease_age <- lapply(plot_dat_list, function(x) { x$dat_cardiovascDisease }) %>%
  dplyr::bind_rows() %>% 
  mutate(diet    = factor(diet,    levels = levels(plot_dat_cardiovascDisease$name)),
         age_cat = factor(age_cat, levels = levels(dat$age_cat)))
plot_dat_weightLoss_age        <- lapply(plot_dat_list, function(x) { x$dat_weightLoss }) %>%
  dplyr::bind_rows() %>% 
  mutate(diet    = factor(diet,    levels = levels(plot_dat_weightLoss$name)),
         age_cat = factor(age_cat, levels = levels(dat$age_cat)))

# create heatmaps
axis_limits <- c(plot_dat_T2Diab_age$share_yes, plot_dat_cardiovascDisease_age$share_yes,
                 plot_dat_weightLoss_age$share_yes) %>% 
  range()

gg1 <- plot_dat_T2Diab_age %>% 
  ggplot(aes(x = age_cat, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  ggtitle("Recommended diet\nfor type-2 diabetes") +
  theme(axis.text.x  = element_text(angle = 45, hjust = 1),
        axis.title.x = element_blank())

gg2 <- plot_dat_cardiovascDisease_age %>% 
  ggplot(aes(x = age_cat, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  ggtitle("Recommended diet\nfor cardiovascular disease") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

gg3 <- plot_dat_weightLoss_age %>% 
  ggplot(aes(x = age_cat, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  ggtitle("Recommended diet\nfor weight loss") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

# joint graphic
(gg1 + gg2 + gg3) +
  patchwork::plot_layout(guides = "collect") +
  patchwork::plot_annotation(title = "Recommended diets by age")

Code
# prepare the datasets
plot_dat_list <- lapply(levels(dat$dietitian_forHowLong_cat), function(x) {
  
  # filter the dataset for dietitians with the specific age category
  dat_practice <- dat %>% filter(dietitian_forHowLong_cat == x)
  
  # T2Diab_diet
  dat_T2Diab <- dat_practice %>% 
    select(starts_with("T2Diab_diet")) %>% 
    select(-T2Diab_diet) %>% 
    pivot_longer(cols = everything(), names_to = "diet") %>% 
    group_by(diet) %>% 
    summarize(share_yes = sum(value == "yes") / n()) %>% 
    arrange(desc(share_yes)) %>% 
    mutate(diet        = gsub(pattern = "T2Diab_diet", replacement = "", x = diet),
           diet        = factor(diet, levels = diet),
           dietitian_forHowLong_cat = x,
           share_label = paste0(round(100 * share_yes), "%")) %>% 
    filter(diet != "NA")

  # cardiovascDisease_diet
  dat_cardiovascDisease <- dat_practice %>% 
    select(starts_with("cardiovascDisease_diet")) %>% 
    select(-cardiovascDisease_diet) %>% 
    pivot_longer(cols = everything(), names_to = "diet") %>% 
    group_by(diet) %>% 
    summarize(share_yes = sum(value == "yes") / n()) %>% 
    arrange(desc(share_yes)) %>% 
    mutate(diet        = gsub(pattern = "cardiovascDisease_diet", replacement = "", x = diet),
           diet        = factor(diet, levels = diet),
           dietitian_forHowLong_cat = x,
           share_label = paste0(round(100 * share_yes), "%")) %>%  
    filter(diet != "NA")
  
  # weightLoss_diet
  dat_weightLoss <- dat_practice %>% 
    select(starts_with("weightLoss_diet")) %>% 
    select(-weightLoss_diet) %>% 
    pivot_longer(cols = everything(), names_to = "diet") %>% 
    group_by(diet) %>% 
    summarize(share_yes = sum(value == "yes") / n()) %>% 
    arrange(desc(share_yes)) %>% 
    mutate(diet        = gsub(pattern = "weightLoss_diet", replacement = "", x = diet),
           diet        = factor(diet, levels = diet),
           dietitian_forHowLong_cat = x,
           share_label = paste0(round(100 * share_yes), "%")) %>%  
    filter(diet != "NA")
  
  list(dat_T2Diab            = dat_T2Diab,
       dat_cardiovascDisease = dat_cardiovascDisease,
       dat_weightLoss        = dat_weightLoss) %>% 
    return()
})

# merge prepared datasets
plot_dat_T2Diab_lon <- lapply(plot_dat_list, function(x) { x$dat_T2Diab }) %>%
  dplyr::bind_rows() %>% 
  mutate(diet                     = factor(diet, levels = levels(plot_dat_T2Diab$name)),
         dietitian_forHowLong_cat = factor(dietitian_forHowLong_cat, levels = levels(dat$dietitian_forHowLong_cat)))
plot_dat_cardiovascDisease_lon <- lapply(plot_dat_list, function(x) { x$dat_cardiovascDisease }) %>%
  dplyr::bind_rows() %>% 
  mutate(diet                     = factor(diet, levels = levels(plot_dat_cardiovascDisease$name)),
         dietitian_forHowLong_cat = factor(dietitian_forHowLong_cat, levels = levels(dat$dietitian_forHowLong_cat)))
plot_dat_weightLoss_lon <- lapply(plot_dat_list, function(x) { x$dat_weightLoss }) %>%
  dplyr::bind_rows() %>% 
  mutate(diet                     = factor(diet, levels = levels(plot_dat_weightLoss$name)),
         dietitian_forHowLong_cat = factor(dietitian_forHowLong_cat, levels = levels(dat$dietitian_forHowLong_cat)))

# create heatmaps
axis_limits <- c(plot_dat_T2Diab_lon$share_yes, plot_dat_cardiovascDisease_lon$share_yes,
                 plot_dat_weightLoss_lon$share_yes) %>% 
  range()

gg1 <- plot_dat_T2Diab_lon %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  ggtitle("Recommended diet\nfor type-2 diabetes") +
  theme(axis.text.x  = element_text(angle = 45, hjust = 1),
        axis.title.x = element_blank())

gg2 <- plot_dat_cardiovascDisease_lon %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  ggtitle("Recommended diet\nfor cardiovascular disease") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

gg3 <- plot_dat_weightLoss_lon %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  ggtitle("Recommended diet\nfor weight loss") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

# joint graphic
(gg1 + gg2 + gg3) +
  patchwork::plot_layout(guides = "collect") +
  patchwork::plot_annotation(title = "Recommended diets by years of practice")

Model-based change

For the WFPB diet, estimate one logistic regression model each for type-2 diabetes, cardiovascular disease and weight loss, on the question ‘yes, I would recommend this diet’ vs. ‘no, I wouldn’t’.

Note: Similar models are not estimated for the vegan diet, as of too small numbers (12 “yes” answers regarding type-2 diabetes, 16 “yes” answers regarding cardiovascular disease, 12 “yes” answers regarding weight loss)

Structure of every logistic regression model:

  • dependent variable: ‘I would recommend this diet for this disease’ (‘yes’ vs. ‘no’)
  • independent variables: subspecialty (one dummy variable for each of the six subspecialties under focus), years of practice (in 10-year-units), education quality on WFPB diets (as numeric score with 0 = ‘No training’, 1 = ‘Very poor’ and 5 = ‘Excellent’)
  • further control variables: education
  • unused control variables: age (since strongly correlated with years of practice), gender (since dietitians are nearly exclusively female), personal diet (since the causality between a dietitian’s personal diet and the response variable is unclear)
Model estimation
Code
model_T2D_WFPB <- estimate_logitModel(dat, y = "T2Diab_dietWFPB")
model_CD_WFPB  <- estimate_logitModel(dat, y = "cardiovascDisease_dietWFPB")
model_WL_WFPB  <- estimate_logitModel(dat, y = "weightLoss_dietWFPB")

AUC type-2 diabetes: 0.56
AUC cardiovascular disease: 0.62
AUC weight loss: 0.61

Results
Code
model_T2D_WFPB$results_tab %>% 
  mutate(exp_estimate = round(exp_estimate, 4),
         exp_CIlower  = round(exp_CIlower, 4),
         exp_CIupper  = round(exp_CIupper, 4)) %>% 
  kable() %>% 
  kable_styling()
param_group parameter exp_estimate exp_CIlower exp_CIupper pvalue
NA (Intercept) 0.2212 0.1026 0.4769 1e-04
area of specialty obesityOrWeightManagementyes 1.6932 0.9051 3.1675 0.0993
area of specialty gastroenterologyyes 1.0747 0.5261 2.1952 0.8433
area of specialty diabetesyes 0.8579 0.4010 1.8353 0.6929
area of specialty paediatricsyes 1.1502 0.5636 2.3474 0.7007
area of specialty elderlyyes 1.2395 0.5904 2.6020 0.5704
area of specialty oncologyyes 0.5917 0.2167 1.6157 0.3059
area of specialty psychEatingyes 2.4544 1.1267 5.3466 0.0238
years of practice dietitian_forHowLong_10yearUnits 0.8998 0.6772 1.1955 0.4663
education educationPostgraduate degree 1.8432 1.0397 3.2674 0.0363
education educationPhD 1.9479 0.5564 6.8193 0.297
WFPB education education_WFPBeducation_num 0.8541 0.6644 1.0980 0.2185
Code
plot_logitRegressionResults(model_T2D_WFPB$results_tab) + ggtitle("Type-2 diabetes")

Code
model_CD_WFPB$results_tab %>% 
  mutate(exp_estimate = round(exp_estimate, 4),
         exp_CIlower  = round(exp_CIlower, 4),
         exp_CIupper  = round(exp_CIupper, 4)) %>% 
  kable() %>% 
  kable_styling()
param_group parameter exp_estimate exp_CIlower exp_CIupper pvalue
NA (Intercept) 0.2951 0.1371 0.6350 0.0018
area of specialty obesityOrWeightManagementyes 2.8326 1.5415 5.2052 8e-04
area of specialty gastroenterologyyes 0.9830 0.4663 2.0720 0.964
area of specialty diabetesyes 0.9963 0.4853 2.0453 0.9919
area of specialty paediatricsyes 0.6147 0.2796 1.3511 0.2259
area of specialty elderlyyes 0.4620 0.1907 1.1194 0.0872
area of specialty oncologyyes 0.4106 0.1361 1.2389 0.1142
area of specialty psychEatingyes 1.9649 0.8659 4.4586 0.1062
years of practice dietitian_forHowLong_10yearUnits 0.8954 0.6743 1.1889 0.445
education educationPostgraduate degree 2.5155 1.4014 4.5151 0.002
education educationPhD 1.2153 0.2941 5.0217 0.7876
WFPB education education_WFPBeducation_num 0.7330 0.5682 0.9456 0.0168
Code
plot_logitRegressionResults(model_CD_WFPB$results_tab)  + ggtitle("Cardiovascular disease")

Code
model_WL_WFPB$results_tab %>% 
  mutate(exp_estimate = round(exp_estimate, 4),
         exp_CIlower  = round(exp_CIlower, 4),
         exp_CIupper  = round(exp_CIupper, 4)) %>% 
  kable() %>% 
  kable_styling()
param_group parameter exp_estimate exp_CIlower exp_CIupper pvalue
NA (Intercept) 0.3524 0.1646 0.7546 0.0073
area of specialty obesityOrWeightManagementyes 2.1367 1.1444 3.9892 0.0172
area of specialty gastroenterologyyes 1.0515 0.5029 2.1985 0.8939
area of specialty diabetesyes 1.0772 0.5152 2.2525 0.8433
area of specialty paediatricsyes 0.9431 0.4460 1.9940 0.8781
area of specialty elderlyyes 0.8495 0.3853 1.8727 0.6858
area of specialty oncologyyes 0.3193 0.0929 1.0976 0.07
area of specialty psychEatingyes 3.3880 1.5412 7.4478 0.0024
years of practice dietitian_forHowLong_10yearUnits 0.7976 0.5931 1.0727 0.1347
education educationPostgraduate degree 1.8127 1.0224 3.2141 0.0418
education educationPhD 0.7675 0.1553 3.7915 0.7454
WFPB education education_WFPBeducation_num 0.7126 0.5507 0.9221 0.01
Code
plot_logitRegressionResults(model_WL_WFPB$results_tab)  + ggtitle("Weight loss")

WFPB: Main nutrients of concern

What is their knowledge of the main nutrients of concern for WFPB diets?

Code
plot_dat_nut <- dat %>% 
  select(starts_with("WFPB_micronutrientsOfConcern_")) %>% 
  pivot_longer(cols = everything()) %>% 
  group_by(name) %>% 
  mutate(sum_yes = sum(value == "yes")) %>% 
  ungroup() %>% 
  arrange(desc(sum_yes)) %>% 
  mutate(name = gsub(pattern = "WFPB_micronutrientsOfConcern_", replacement = "", name),
         name = factor(name, levels = rev(unique(name))))
plot_dat_nut %>% 
  ggplot(aes(y = name, fill = value)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("lightgray", unname(col_vector["knowledge"]))) +
  scale_x_continuous(labels = scales::label_percent()) +
  xlab("Agreement [%]") +
  ggtitle("Micronutrients of concern") +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

Univariate description

by subspecialty, age, years of practice.

Code
# prepare the dataset
plot_dat_list <- lapply(1:length(subspecialties_varVector), function(i) {
  
  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # filter the dataset for dietitians with the specific subspecialty
  dat_subsp <- dat
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  results_dat <- dat_subsp %>% 
    select(starts_with("WFPB_micronutrientsOfConcern_")) %>% 
    pivot_longer(cols = everything(), names_to = "micronutrient") %>% 
    group_by(micronutrient) %>% 
    summarize(share_yes = sum(value == "yes") / n()) %>% 
    arrange(desc(share_yes)) %>% 
    mutate(micronutrient = gsub(pattern = "WFPB_micronutrientsOfConcern_", replacement = "", x = micronutrient),
           micronutrient = factor(micronutrient, levels = rev(micronutrient)),
           subspecialty  = subsp_label,
           share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat <- plot_dat_list %>% bind_rows() %>% 
  mutate(micronutrient = factor(micronutrient, levels = levels(plot_dat_nut$name)),
         subspecialty  = factor(subspecialty,  levels = names(subspecialties_varVector)))

# heatmap
plot_dat %>% 
  ggplot(aes(x = subspecialty, y = micronutrient, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      labels = scales::label_percent()) +
  ggtitle("WFPB micronutrients of concern\nby subspecialties") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Code
# prepare the dataset
plot_dat_list <- lapply(levels(dat$age_cat), function(x) {
  
  # filter the dataset for dietitians with the specific age group
  dat_age <- dat %>% filter(age_cat == x)
  
  results_dat <- dat_age %>% 
    select(starts_with("WFPB_micronutrientsOfConcern_")) %>% 
    pivot_longer(cols = everything(), names_to = "micronutrient") %>% 
    group_by(micronutrient) %>% 
    summarize(share_yes = sum(value == "yes") / n()) %>% 
    arrange(desc(share_yes)) %>% 
    mutate(micronutrient = gsub(pattern = "WFPB_micronutrientsOfConcern_", replacement = "", x = micronutrient),
           micronutrient = factor(micronutrient, levels = rev(micronutrient)),
           age_cat       = x,
           share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat <- plot_dat_list %>% bind_rows() %>% 
  mutate(micronutrient = factor(micronutrient, levels = levels(plot_dat_nut$name)),
         age_cat       = factor(age_cat,       levels = levels(dat$age_cat))) 

# heatmap
plot_dat %>% 
  ggplot(aes(x = age_cat, y = micronutrient, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      labels = scales::label_percent()) +
  ggtitle("WFPB micronutrients of concern\nby age groups") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Code
# prepare the dataset
plot_dat_list <- lapply(levels(dat$dietitian_forHowLong_cat), function(x) {
  
  # filter the dataset for dietitians with the specific age group
  dat_practice <- dat %>% filter(dietitian_forHowLong_cat == x)
  
  results_dat <- dat_practice %>% 
    select(starts_with("WFPB_micronutrientsOfConcern_")) %>% 
    pivot_longer(cols = everything(), names_to = "micronutrient") %>% 
    group_by(micronutrient) %>% 
    summarize(share_yes = sum(value == "yes") / n()) %>% 
    arrange(desc(share_yes)) %>% 
    mutate(micronutrient = gsub(pattern = "WFPB_micronutrientsOfConcern_", replacement = "", x = micronutrient),
           micronutrient = factor(micronutrient, levels = rev(micronutrient)),
           dietitian_forHowLong_cat = x,
           share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat <- plot_dat_list %>% bind_rows() %>% 
  mutate(micronutrient            = factor(micronutrient, levels = levels(plot_dat_nut$name)),
         dietitian_forHowLong_cat = factor(dietitian_forHowLong_cat, levels = levels(dat$dietitian_forHowLong_cat))) 

# heatmap
plot_dat %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = micronutrient, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      labels = scales::label_percent()) +
  ggtitle("WFPB micronutrients of concern\nby dietitians' years of practice") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Model-based change

Individually for every nutrient apart from potassium (which has too few “yes” answers): Logistic regression on the question ‘yes, that’s a nutrient of concern’ vs. ‘no, it isn’t’

Structure of every logistic regression model:

  • dependent variable: ‘This is a WFPB micronutrient is of concern’ (‘yes’ vs. ‘no’)
  • independent variables: subspecialty (one dummy variable for each of the six subspecialties under focus), years of practice (in 10-year-units), education quality on WFPB diets (as numeric score with 0 = ‘No training’, 1 = ‘Very poor’ and 5 = ‘Excellent’)
  • further control variables: education
  • unused control variables: age (since strongly correlated with years of practice), gender (since dietitians are nearly exclusively female), personal diet (since the causality between a dietitian’s personal diet and the response variable is unclear)
Model estimation
Code
model_vitB12  <- estimate_logitModel(dat, y = "WFPB_micronutrientsOfConcern_vitB12")
model_iron    <- estimate_logitModel(dat, y = "WFPB_micronutrientsOfConcern_iron")
model_calcium <- estimate_logitModel(dat, y = "WFPB_micronutrientsOfConcern_calcium")
model_iodine  <- estimate_logitModel(dat, y = "WFPB_micronutrientsOfConcern_iodine")
model_LComg3  <- estimate_logitModel(dat, y = "WFPB_micronutrientsOfConcern_LComega3")
model_vitD    <- estimate_logitModel(dat, y = "WFPB_micronutrientsOfConcern_vitD")
model_zinc    <- estimate_logitModel(dat, y = "WFPB_micronutrientsOfConcern_zinc")
model_SComg3  <- estimate_logitModel(dat, y = "WFPB_micronutrientsOfConcern_SComega3")
model_folate  <- estimate_logitModel(dat, y = "WFPB_micronutrientsOfConcern_folate")
model_selen   <- estimate_logitModel(dat, y = "WFPB_micronutrientsOfConcern_selenium")
model_thiam   <- estimate_logitModel(dat, y = "WFPB_micronutrientsOfConcern_thiamine")
model_choline <- estimate_logitModel(dat, y = "WFPB_micronutrientsOfConcern_choline")

AUC vitamin B12: 0.61
AUC iron: 0.47
AUC calcium: 0.53
AUC iodine: 0.52
AUC LComega3: 0.52
AUC vitamin D: 0.58
AUC zinc: 0.55
AUC SComega3: 0.59
AUC folate: 0.33
AUC selenium: 0.36
AUC thiamine: 0.34
AUC choline: 0.3

Effect plots
Code
plot_logitRegressionResults(model_vitB12$results_tab) + ggtitle("Vitamin B12")

Code
plot_logitRegressionResults(model_iron$results_tab) + ggtitle("Iron")

Code
plot_logitRegressionResults(model_calcium$results_tab) + ggtitle("Calcium")

Code
plot_logitRegressionResults(model_iodine$results_tab) + ggtitle("Iodine")

Code
plot_logitRegressionResults(model_LComg3$results_tab) + ggtitle("LComega3")

Code
plot_logitRegressionResults(model_vitD$results_tab) + ggtitle("Vitamin D")

Code
plot_logitRegressionResults(model_zinc$results_tab) + ggtitle("Zinc")

Code
plot_logitRegressionResults(model_SComg3$results_tab) + ggtitle("SComega3")

Code
plot_logitRegressionResults(model_folate$results_tab) + ggtitle("Folate")
Warning in scale_y_continuous("Odds Ratio, on log2-transformed scale", trans =
"log2"): log-2 transformation introduced infinite values.

Code
plot_logitRegressionResults(model_selen$results_tab) + ggtitle("Selenium")

Code
plot_logitRegressionResults(model_thiam$results_tab) + ggtitle("Thiamine")
Warning in scale_y_continuous("Odds Ratio, on log2-transformed scale", trans =
"log2"): log-2 transformation introduced infinite values.

Code
plot_logitRegressionResults(model_choline$results_tab) + ggtitle("Choline")

WFPB: Increased risk of malnutrition or eating disorders

Do they think WFPB diets increase the risk of of i) malnutrition or ii) eating disorders?

Code
dat %>% 
  select(Pbdiet_increasedRiskForMalnutrition, Pbdiet_increasedRiskEatingDisorder) %>% 
  pivot_longer(cols = everything(), names_to = "item") %>% 
  ggplot(aes(y = item, fill = value)) +
  geom_bar(position = "fill") +
  scale_x_continuous("Relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "PRGn") +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

Plant proteins: incomplete

Do they think plant-proteins are incomplete?

Code
dat %>% 
  ggplot(aes(y = "", fill = plantProtein_isIncomplete)) +
  geom_bar(position = "fill") +
  scale_x_continuous("Relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "PRGn", direction = -1) +
  ggtitle("Plant proteins are an incomplete source of protein") +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

Univariate description

by subspecialty, age, years of practice

Code
# prepare the dataset
plot_dat_list <- lapply(1:length(subspecialties_varVector), function(i) {
  
  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # filter the dataset for dietitians with the specific subspecialty
  dat_subsp <- dat
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  results_dat <- data.frame(plantProtein_isIncomplete = levels(dat_subsp$plantProtein_isIncomplete),
             share_yes                   = prop.table(table(dat_subsp$plantProtein_isIncomplete)) %>% as.vector,
             subspecialty                = subsp_label) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat_proteinSubsp <- plot_dat_list %>% bind_rows() %>% 
  mutate(subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)),
         plantProtein_isIncomplete = factor(plantProtein_isIncomplete, levels = levels(dat$plantProtein_isIncomplete)))

# heatmap
plot_dat_proteinSubsp %>% 
  ggplot(aes(x = subspecialty, y = plantProtein_isIncomplete, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      labels = scales::label_percent()) +
  ggtitle("Plant proteins are an incomplete source of protein\nby subspecialties") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Code
# prepare the dataset
plot_dat_list <- lapply(levels(dat$age_cat), function(x) {
  
  # filter the dataset for dietitians with the specific age group
  dat_age <- dat %>% filter(age_cat == x)
  
  results_dat <- data.frame(plantProtein_isIncomplete = levels(dat_age$plantProtein_isIncomplete),
                            share_yes                   = prop.table(table(dat_age$plantProtein_isIncomplete)) %>% as.vector,
                            age_cat                     = x) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat_proteinAge <- plot_dat_list %>% bind_rows() %>% 
  mutate(age_cat                     = factor(age_cat, levels = levels(dat$age_cat)),
         plantProtein_isIncomplete  = factor(plantProtein_isIncomplete, levels = levels(dat$plantProtein_isIncomplete)))

# heatmap
plot_dat_proteinAge %>% 
  ggplot(aes(x = age_cat, y = plantProtein_isIncomplete, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      labels = scales::label_percent()) +
  ggtitle("Plant proteins are an incomplete source of protein\nby age groups") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Code
# prepare the dataset
plot_dat_list <- lapply(levels(dat$dietitian_forHowLong_cat), function(x) {
  
  # filter the dataset for dietitians with the specific years of practice
  dat_practice <- dat %>% filter(dietitian_forHowLong_cat == x)
  
  results_dat <- data.frame(plantProtein_isIncomplete  = levels(dat_practice$plantProtein_isIncomplete),
                            share_yes                   = prop.table(table(dat_practice$plantProtein_isIncomplete)) %>% as.vector,
                            dietitian_forHowLong_cat    = x) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat_proteinPractice <- plot_dat_list %>% bind_rows() %>% 
  mutate(dietitian_forHowLong_cat    = factor(dietitian_forHowLong_cat, levels = levels(dat$dietitian_forHowLong_cat)),
         plantProtein_isIncomplete  = factor(plantProtein_isIncomplete, levels = levels(dat$plantProtein_isIncomplete)))

# heatmap
plot_dat_proteinPractice %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = plantProtein_isIncomplete, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      labels = scales::label_percent()) +
  ggtitle("Plant proteins are an incomplete source of protein\nby dietitians' years of practice") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Model-based change

Estimation of a Logistic regression model:

  • dependent variable: binarized incomplete plant protein item (‘Yes, (strongly) agree’ vs. ‘Not sure / No, (strongly) disagree’)
  • independent variables: subspecialty (one dummy variable for each of the six subspecialties under focus), years of practice (in 10-year-units), education quality on WFPB diets (as numeric score with 0 = ‘No training’, 1 = ‘Very poor’ and 5 = ‘Excellent’)
  • further control variables: education
  • unused control variables: age (since strongly correlated with years of practice), gender (since dietitians are nearly exclusively female), personal diet (since the causality between a dietitian’s personal diet and the response variable is unclear)
Code
# calculate the response variable
dat <- dat %>% 
  mutate(plantProteinIsIncomplete_binary = case_when(grepl("Yes", plantProtein_isIncomplete) ~ "Yes, (strongly) agree",
                                                     plantProtein_isIncomplete == "Not sure" ~ "Not sure / No, (strongly) disagree",
                                                     grepl("No", plantProtein_isIncomplete)  ~ "Not sure / No, (strongly) disagree",
                                                     TRUE                                    ~ NA_character_),
         plantProteinIsIncomplete_binary = factor(plantProteinIsIncomplete_binary, levels = c("Not sure / No, (strongly) disagree","Yes, (strongly) agree")))

# model estimation
model <- gam(plantProteinIsIncomplete_binary ~ 
               dieteticsSubspecialty_obesityOrWeightManagement +
               dieteticsSubspecialty_gastroenterology +
               dieteticsSubspecialty_diabetes +
               dieteticsSubspecialty_paediatrics +
               dieteticsSubspecialty_elderly +
               dieteticsSubspecialty_oncology +
               dieteticsSubspecialty_psychEating +
               dietitian_forHowLong_10yearUnits +
               education +
               education_WFPBeducation_num,
             family = binomial(link = "logit"),
             data = dat)

# prepare the results table
results_tab_plantProtein <- summary(model)$p.table %>% 
  as.data.frame() %>% 
  mutate(parameter = row.names(.)) %>% 
  rename(estimate = "Estimate",
         se       = "Std. Error",
         pvalue   = "Pr(>|z|)") %>% 
  mutate(CI_lower     = estimate - qnorm(.975) * se,
         CI_upper     = estimate + qnorm(.975) * se,
         exp_estimate = exp(estimate),
         exp_CIlower  = exp(CI_lower),
         exp_CIupper  = exp(CI_upper),
         pvalue       = case_when(pvalue < .0001 ~ "<.0001",
                                  TRUE           ~ as.character(round(pvalue, 4))),
         param_group  = case_when(grepl("Subspecialty", parameter)  ~ "area of specialty",
                                  grepl("forHowLong", parameter)    ~ "years of practice",
                                  grepl("WFPBeducation", parameter) ~ "WFPB education",
                                  grepl("education", parameter)     ~ "education"),
         param_group  = factor(param_group, levels = c("area of specialty", "years of practice", "education", "WFPB education")),
         parameter    = gsub(parameter, pattern = "dieteticsSubspecialty_", replacement = "")) %>% 
  select(param_group, parameter, exp_estimate, exp_CIlower, exp_CIupper, pvalue)
row.names(results_tab_plantProtein) <- NULL
results_tab_plantProtein %>% 
  mutate(exp_estimate = round(exp_estimate, 4),
         exp_CIlower  = round(exp_CIlower, 4),
         exp_CIupper  = round(exp_CIupper, 4)) %>% 
  kable() %>% 
  kable_styling()
param_group parameter exp_estimate exp_CIlower exp_CIupper pvalue
NA (Intercept) 2.3644 1.1596 4.8211 0.0179
area of specialty obesityOrWeightManagementyes 0.6626 0.3599 1.2197 0.1862
area of specialty gastroenterologyyes 2.4942 1.1023 5.6436 0.0282
area of specialty diabetesyes 0.8775 0.4381 1.7574 0.7122
area of specialty paediatricsyes 1.4394 0.6873 3.0141 0.3342
area of specialty elderlyyes 0.9675 0.4606 2.0323 0.9306
area of specialty oncologyyes 0.4275 0.1980 0.9230 0.0305
area of specialty psychEatingyes 1.4424 0.5600 3.7154 0.448
years of practice dietitian_forHowLong_10yearUnits 0.9041 0.6960 1.1744 0.4499
education educationPostgraduate degree 0.8589 0.5015 1.4709 0.5795
education educationPhD 0.5294 0.1704 1.6452 0.2716
WFPB education education_WFPBeducation_num 1.3335 1.0424 1.7060 0.022
Code
# effect plot
results_tab_plantProtein %>% 
  filter(parameter != "(Intercept)") %>% 
  ggplot(aes(x = parameter)) +
  geom_hline(yintercept = 1, lty = 2, col = "gray50") +
  geom_pointrange(aes(y = exp_estimate, ymin = exp_CIlower, ymax = exp_CIupper, col = param_group)) +
  scale_y_continuous("Odds Ratio, on log2-transformed scale", trans = "log2") +
  facet_wrap(~ param_group, scales = "free_x", nrow = 1) +
  theme(axis.text.x     = element_text(angle = 45, hjust = 1),
        axis.title.x    = element_blank(),
        legend.position = "none")

The model’s goodness-of-fit is determined by calculating the AUC (Area Under the Curve). Therefore, the dataset is first randomly split into a training dataset (80% of the data) and a test dataset (20% of the data). The regression model is then re-estimated on the training data, and the AUC is calculated on the test data.

Code
# basis: split data into a training set (80% of the data) and a test set (20%)
set.seed(2024)
train_obs <- sample(1:nrow(dat), size = round(0.8*nrow(dat)), replace = FALSE)
dat_train <- dat %>% slice(train_obs)
dat_test  <- dat %>% slice(-train_obs)
Code
# 1) estimate the model on the training set
model_train <- gam(plantProteinIsIncomplete_binary ~ 
                     dieteticsSubspecialty_obesityOrWeightManagement +
                     dieteticsSubspecialty_gastroenterology +
                     dieteticsSubspecialty_diabetes +
                     dieteticsSubspecialty_paediatrics +
                     dieteticsSubspecialty_elderly +
                     dieteticsSubspecialty_oncology +
                     dieteticsSubspecialty_psychEating +
                     dietitian_forHowLong_10yearUnits +
                     education,
                   family = binomial(link = "logit"),
                   data = dat_train)
# 2) calculate ROC curve and the corresponding AUC value
pred <- data.frame("outcome" = as.integer(dat_test$plantProteinIsIncomplete_binary) - 1,
                   "pred"    = predict(model_train, newdata = dat_test))
gg <- ggplot(pred, aes(m = pred, d = outcome)) + geom_roc() + ggtitle("ROC curve")
auc <- calc_auc(gg)$AUC

AUC: 0.64

WFPB: Life-cycle knowledge score

Does knowledge of WFPB diets in the lifecycle change or is determined by i) age or ii) years practising as a dietitian?

We define a life-cycle WFPB knowledge score based on these five items:

  1. A well-planned whole food plant-based diet is suitable and healthy for all life stages.
  2. A well-planned whole food plant-based diet is suitable in all stages of pregnancy and lactation.
  3. It is possible for children (infants and toddlers) to meet all nutritional requirements on a well-planned whole food plant-based diet?
  4. It is possible for adolescents and teenagers to meet all nutritional requirements on a well-planned whole food 5. plant-based diet?
  5. It is difficult for older persons to achieve their energy and protein requirements on a well-planned, whole food plant-based diet?

Based on these five items, the knowledge score for every person is calculated by taking the average of a person’s responses to these items (with the lowest response category ‘No, strongly disagree’ encoded as 0, and the highest ‘Yes, strongly agree’ encoded as 1). Categories ‘Disagree’, ‘Agree’ and ‘Unsure’ are set to values ‘0’ (since it’s a wrong answer), ‘0.75’, and ‘0.75 / 2 = 0.375’ (to reflect the middle between agree and disagree), respectively.
Since item number five is encoded in different direction (agreement = negative view of WFPB diets) compared to items one to four (agreement = positive view of WFPB diets), the direction of the Likert scale of item number five is turned around before calculating the knowledge score.

Code
dat <- dat %>% 
  mutate(item1_num = case_when(grepl("disagree", wpWFPB_isHealthy_allLifeStages)       ~ 0,
                               wpWFPB_isHealthy_allLifeStages == "Not sure"            ~ 0.375,
                               wpWFPB_isHealthy_allLifeStages == "Yes, agree"          ~ 0.75,
                               wpWFPB_isHealthy_allLifeStages == "Yes, strongly agree" ~ 1),
         item2_num = case_when(grepl("disagree", wpWFPB_isHealthy_pregnancyLact)       ~ 0,
                               wpWFPB_isHealthy_pregnancyLact == "Not sure"            ~ 0.375,
                               wpWFPB_isHealthy_pregnancyLact == "Yes, agree"          ~ 0.75,
                               wpWFPB_isHealthy_pregnancyLact == "Yes, strongly agree" ~ 1),
         item3_num = case_when(grepl("disagree", wpWFPB_isHealthy_children)       ~ 0,
                               wpWFPB_isHealthy_children == "Not sure"            ~ 0.375,
                               wpWFPB_isHealthy_children == "Yes, agree"          ~ 0.75,
                               wpWFPB_isHealthy_children == "Yes, strongly agree" ~ 1),
         item4_num = case_when(grepl("disagree", wpWFPB_isHealthy_teens)       ~ 0,
                               wpWFPB_isHealthy_teens == "Not sure"            ~ 0.375,
                               wpWFPB_isHealthy_teens == "Yes, agree"          ~ 0.75,
                               wpWFPB_isHealthy_teens == "Yes, strongly agree" ~ 1),
         item5_num = case_when(wpWFPB_isDifficult_older == "No, strongly disagree" ~ 1,
                               wpWFPB_isDifficult_older == "No, disagree"          ~ 0.75,
                               wpWFPB_isDifficult_older == "Not sure"              ~ 0.375,
                               grepl("agree", wpWFPB_isDifficult_older)            ~ 0)) %>% 
  rowwise() %>% 
  mutate(knowledge_score_lifeCycle = mean(c(item1_num, item2_num, item3_num, item4_num, item5_num))) %>% 
  ungroup()

# boxplot
ggplot(dat, aes(x = "", y = knowledge_score_lifeCycle)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Life-cycle knowledge score") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Code
# alternative barplot
ggplot(dat, aes(x = knowledge_score_lifeCycle, y = after_stat(count / sum(count)))) +
  geom_bar(fill = col_vector["knowledge"]) +
  scale_y_continuous("relative frequency", labels = scales::label_percent()) +
  scale_x_continuous("life cycle knowledge score", labels = scales::label_percent(),
                     limits = c(0,NA)) +
  theme(panel.grid.minor = element_blank())

Univariate description

Some general statistics:

  • arithmetic mean score (min – max): 61% (0% – 100%)
  • quantiles: 0% (2.5% quantile), 15% (5% quantile), 45% (25% quantile), 62% (50% quantile), 80% (75% quantile), 100% (95% quantile), 100% (97.5% quantile)
  • 7% of the respondents scored 100%
  • 37% of the respondents scored 75% or higher
  • 71% of the respondents scored 50% or higher
  • 91% of the respondents scored 25% or higher
  • 3% of the respondents scored 0%

The following plots are conditioned on subspecialty, age, years of practice

Code
plotDat_list <- lapply(1:length(subspecialties_varVector), function(i) {

  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # rename the respective subspecialty column, for easier handling
  dat_subsp <- dat
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  plot_dat <- data.frame(subspecialty              = subsp_label,
                         knowledge_score_lifeCycle = dat_subsp$knowledge_score_lifeCycle)
  
  return(plot_dat)
})
plot_dat <- plotDat_list %>% bind_rows() %>% 
  mutate(subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)))

ggplot(plot_dat, aes(x = subspecialty, y = knowledge_score_lifeCycle)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Life-cycle knowledge score\nby subspecialty") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Code
dat %>% 
  filter(!is.na(age_cat)) %>% 
  ggplot(aes(x = age_cat, y = knowledge_score_lifeCycle)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Life-cycle knowledge score\nby age") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Code
dat %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = knowledge_score_lifeCycle)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Life-cycle knowledge score\nby years of practice") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Model-based change

Estimation of a Gaussian regression model:

  • dependent variable: Life-cycle knowledge score (quasi-metric, assumed as metric here)
  • independent variables: subspecialty (one dummy variable for each of the six subspecialties under focus), years of practice (in 10-year-units), education quality on WFPB diets (as numeric score with 0 = ‘No training’, 1 = ‘Very poor’ and 5 = ‘Excellent’)
  • further control variables: education
  • unused control variables: age (since strongly correlated with years of practice), gender (since dietitians are nearly exclusively female), personal diet (since the causality between a dietitian’s personal diet and the response variable is unclear)
Code
# model estimation
model <- gam(knowledge_score_lifeCycle ~ 
               dieteticsSubspecialty_obesityOrWeightManagement +
               dieteticsSubspecialty_gastroenterology +
               dieteticsSubspecialty_diabetes +
               dieteticsSubspecialty_paediatrics +
               dieteticsSubspecialty_elderly +
               dieteticsSubspecialty_oncology +
               dieteticsSubspecialty_psychEating +
               dietitian_forHowLong_10yearUnits +
               education +
               education_WFPBeducation_num,
             data = dat)

# prepare the results table
results_tab_LCscore <- summary(model)$p.table %>% 
  as.data.frame() %>% 
  mutate(parameter = row.names(.)) %>% 
  rename(estimate = "Estimate",
         se       = "Std. Error",
         pvalue   = "Pr(>|t|)") %>% 
  mutate(CI_lower = estimate - qnorm(.975) * se,
         CI_upper = estimate + qnorm(.975) * se,
         pvalue   = case_when(pvalue < .0001 ~ "<.0001",
                              TRUE           ~ as.character(round(pvalue, 4))),
         param_group = case_when(grepl("Subspecialty", parameter)  ~ "area of specialty",
                                 grepl("forHowLong", parameter)    ~ "years of practice",
                                 grepl("WFPBeducation", parameter) ~ "WFPB education",
                                 grepl("education", parameter)     ~ "education"),
         param_group = factor(param_group, levels = c("area of specialty", "years of practice", "education", "WFPB education")),
         parameter    = gsub(parameter, pattern = "dieteticsSubspecialty_", replacement = "")) %>% 
  select(param_group, parameter, estimate, CI_lower, CI_upper, pvalue)
row.names(results_tab_LCscore) <- NULL
results_tab_LCscore %>% 
  mutate(estimate = round(estimate, 4),
         CI_lower = round(CI_lower, 4),
         CI_upper = round(CI_upper, 4)) %>% 
  kable() %>% 
  kable_styling()
param_group parameter estimate CI_lower CI_upper pvalue
NA (Intercept) 0.6453 0.5682 0.7224 <.0001
area of specialty obesityOrWeightManagementyes 0.0511 -0.0160 0.1182 0.1366
area of specialty gastroenterologyyes -0.0311 -0.1038 0.0416 0.403
area of specialty diabetesyes 0.0004 -0.0772 0.0780 0.9915
area of specialty paediatricsyes 0.0032 -0.0708 0.0772 0.9326
area of specialty elderlyyes -0.0061 -0.0833 0.0711 0.8772
area of specialty oncologyyes -0.0888 -0.1773 -0.0003 0.05
area of specialty psychEatingyes 0.0578 -0.0317 0.1472 0.2064
years of practice dietitian_forHowLong_10yearUnits -0.0174 -0.0457 0.0108 0.2277
education educationPostgraduate degree 0.0037 -0.0528 0.0602 0.8984
education educationPhD -0.0680 -0.1991 0.0632 0.3104
WFPB education education_WFPBeducation_num -0.0094 -0.0350 0.0163 0.4763
Code
# effect plot
results_tab_LCscore %>% 
  filter(parameter != "(Intercept)") %>% 
  ggplot(aes(x = parameter)) +
  geom_hline(yintercept = 0, lty = 2, col = "gray50") +
  geom_pointrange(aes(y = estimate, ymin = CI_lower, ymax = CI_upper, col = param_group)) +
  facet_wrap(~ param_group, scales = "free_x", nrow = 1) +
  theme(axis.text.x     = element_text(angle = 45, hjust = 1),
        axis.title.x    = element_blank(),
        legend.position = "none")

R²: 0.004
Share of explained deviance: 3.8%

Residual structure looks acceptable.

Code
# # check the residuals
# resid_dat <- data.frame(prediction = model$fitted.values,
#                         residuum   = model$residuals)
# ggplot(resid_dat, aes(x = prediction, y = residuum)) +
#   geom_point() +
#   geom_smooth(se = FALSE) +
#   ggtitle("Model residuals vs. predicted values")

WFPB: Critical Micronutrients knowledge score

Does knowledge about critical micronutrients in WFPB diets change or is determined by i) age or ii) years practising as a dietitian?

We define a critical micronutrient WFPB knowledge score based on these seven micronutrients:

  1. Vitamin B12
  2. Calcium
  3. Vitamin D
  4. Zinc
  5. Long-chain omega 3
  6. Iodine
  7. Selenium

Based on these seven micronutrients, the knowledge score for every person is calculated by taking the average of a person’s responses to these items. The obtained score then lies between zero and one and reflects the share of the above micronutrients which were correctly answered as being critical.

Code
dat <- dat %>% 
  mutate(item1_num = as.numeric(WFPB_micronutrientsOfConcern_vitB12) - 1,
         item2_num = as.numeric(WFPB_micronutrientsOfConcern_calcium) - 1,
         item3_num = as.numeric(WFPB_micronutrientsOfConcern_vitD) - 1,
         item4_num = as.numeric(WFPB_micronutrientsOfConcern_zinc) - 1,
         item5_num = as.numeric(WFPB_micronutrientsOfConcern_LComega3) - 1,
         item6_num = as.numeric(WFPB_micronutrientsOfConcern_iodine) - 1,
         item7_num = as.numeric(WFPB_micronutrientsOfConcern_selenium) - 1) %>% 
  rowwise() %>% 
  mutate(knowledge_score_criticalNutrients = mean(c(item1_num, item2_num, item3_num, item4_num, item5_num, item6_num, item7_num))) %>% 
  ungroup()

# boxplot
ggplot(dat, aes(x = "", y = knowledge_score_criticalNutrients)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Critical micronutrients knowledge score") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Code
# alternative barplot
ggplot(dat, aes(x = knowledge_score_criticalNutrients, y = after_stat(count / sum(count)))) +
  geom_bar(fill = col_vector["knowledge"]) +
  scale_y_continuous("relative frequency", labels = scales::label_percent()) +
  scale_x_continuous("critical micronutrients knowledge score",
                     breaks = (0:7)/7, labels = paste0(0:7, "/7")) +
  theme(panel.grid.minor = element_blank())

Univariate description

Some general statistics:

  • arithmetic mean score (min – max): 50% (0% – 100%)
  • quantiles: 14% (2.5% quantile), 14% (5% quantile), 29% (25% quantile), 43% (50% quantile), 71% (75% quantile), 86% (95% quantile), 86% (97.5% quantile)
  • 2% of the respondents scored 100%
  • 12% of the respondents scored 75% or higher
  • 48% of the respondents scored 50% or higher
  • 89% of the respondents scored 25% or higher
  • 1% of the respondents scored 0%

The following plots are conditioned on subspecialty, age, years of practice

Code
plotDat_list <- lapply(1:length(subspecialties_varVector), function(i) {

  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # rename the respective subspecialty column, for easier handling
  dat_subsp <- dat
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  plot_dat <- data.frame(subspecialty              = subsp_label,
                         knowledge_score_criticalNutrients = dat_subsp$knowledge_score_criticalNutrients)
  
  return(plot_dat)
})
plot_dat <- plotDat_list %>% bind_rows() %>% 
  mutate(subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)))

ggplot(plot_dat, aes(x = subspecialty, y = knowledge_score_criticalNutrients)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Critical micronutrients knowledge score\nby subspecialty") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Code
dat %>% 
  filter(!is.na(age_cat)) %>% 
  ggplot(aes(x = age_cat, y = knowledge_score_criticalNutrients)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Critical micronutrients knowledge score\nby age") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Code
dat %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = knowledge_score_criticalNutrients)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Critical micronutrients knowledge score\nby years of practice") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Model-based change

Estimation of a Gaussian regression model:

  • dependent variable: Critical micronutrients knowledge score (quasi-metric, assumed as metric here)
  • independent variables: subspecialty (one dummy variable for each of the six subspecialties under focus), years of practice (in 10-year-units), education quality on WFPB diets (as numeric score with 0 = ‘No training’, 1 = ‘Very poor’ and 5 = ‘Excellent’)
  • further control variables: education
  • unused control variables: age (since strongly correlated with years of practice), gender (since dietitians are nearly exclusively female), personal diet (since the causality between a dietitian’s personal diet and the response variable is unclear)
Code
# model estimation
model <- gam(knowledge_score_criticalNutrients ~ 
               dieteticsSubspecialty_obesityOrWeightManagement +
               dieteticsSubspecialty_gastroenterology +
               dieteticsSubspecialty_diabetes +
               dieteticsSubspecialty_paediatrics +
               dieteticsSubspecialty_elderly +
               dieteticsSubspecialty_oncology +
               dieteticsSubspecialty_psychEating +
               dietitian_forHowLong_10yearUnits +
               education +
               education_WFPBeducation_num,
             data = dat)

# prepare the results table
results_tab_CNscore <- summary(model)$p.table %>% 
  as.data.frame() %>% 
  mutate(parameter = row.names(.)) %>% 
  rename(estimate = "Estimate",
         se       = "Std. Error",
         pvalue   = "Pr(>|t|)") %>% 
  mutate(CI_lower = estimate - qnorm(.975) * se,
         CI_upper = estimate + qnorm(.975) * se,
         pvalue   = case_when(pvalue < .0001 ~ "<.0001",
                              TRUE           ~ as.character(round(pvalue, 4))),
         param_group = case_when(grepl("Subspecialty", parameter)  ~ "area of specialty",
                                 grepl("forHowLong", parameter)    ~ "years of practice",
                                 grepl("WFPBeducation", parameter) ~ "WFPB education",
                                 grepl("education", parameter)     ~ "education"),
         param_group = factor(param_group, levels = c("area of specialty", "years of practice", "education", "WFPB education")),
         parameter    = gsub(parameter, pattern = "dieteticsSubspecialty_", replacement = "")) %>% 
  select(param_group, parameter, estimate, CI_lower, CI_upper, pvalue)
row.names(results_tab_CNscore) <- NULL
results_tab_CNscore %>% 
  mutate(estimate = round(estimate, 4),
         CI_lower = round(CI_lower, 4),
         CI_upper = round(CI_upper, 4)) %>% 
  kable() %>% 
  kable_styling()
param_group parameter estimate CI_lower CI_upper pvalue
NA (Intercept) 0.4913 0.4222 0.5603 <.0001
area of specialty obesityOrWeightManagementyes 0.0046 -0.0555 0.0647 0.8811
area of specialty gastroenterologyyes 0.0166 -0.0485 0.0818 0.6167
area of specialty diabetesyes -0.0080 -0.0775 0.0614 0.8206
area of specialty paediatricsyes 0.0693 0.0031 0.1356 0.0411
area of specialty elderlyyes 0.0181 -0.0510 0.0872 0.6083
area of specialty oncologyyes -0.0734 -0.1527 0.0058 0.0704
area of specialty psychEatingyes 0.0462 -0.0339 0.1263 0.2589
years of practice dietitian_forHowLong_10yearUnits 0.0010 -0.0243 0.0263 0.937
education educationPostgraduate degree 0.0089 -0.0417 0.0595 0.7304
education educationPhD 0.0318 -0.0856 0.1493 0.5958
WFPB education education_WFPBeducation_num -0.0074 -0.0304 0.0156 0.5289
Code
# effect plot
results_tab_CNscore %>% 
  filter(parameter != "(Intercept)") %>% 
  ggplot(aes(x = parameter)) +
  geom_hline(yintercept = 0, lty = 2, col = "gray50") +
  geom_pointrange(aes(y = estimate, ymin = CI_lower, ymax = CI_upper, col = param_group)) +
  facet_wrap(~ param_group, scales = "free_x", nrow = 1) +
  theme(axis.text.x     = element_text(angle = 45, hjust = 1),
        axis.title.x    = element_blank(),
        legend.position = "none")

R²: -0.001
Share of explained deviance: 3.3%

Residual structure looks acceptable.

Code
# # check the residuals
# resid_dat <- data.frame(prediction = model$fitted.values,
#                         residuum   = model$residuals)
# ggplot(resid_dat, aes(x = prediction, y = residuum)) +
#   geom_point() +
#   geom_smooth(se = FALSE) +
#   ggtitle("Model residuals vs. predicted values")

WFPB: Uncritical Micronutrients knowledge score

Does knowledge about uncritical micronutrients in WFPB diets change or is determined by i) age or ii) years practising as a dietitian?

We define an uncritical micronutrient WFPB knowledge score based on these six micronutrients:

  1. Iron
  2. Folic acid
  3. Choline
  4. Thiamine
  5. Potassium
  6. Short-chain omega 3

Based on these six micronutrients, the knowledge score for every person is calculated by taking the average of a person’s responses to these items. This score is then inverted (i.e., ‘new_score = 1 - old_score’), s.t. value 1 reflects all correct answers, and value 0 all wrong answers. The obtained score lies between zero and one and reflects the share of the above micronutrients which were correctly answered as not being critical.

Code
dat <- dat %>% 
  mutate(item1_num = as.numeric(WFPB_micronutrientsOfConcern_iron) - 1,
         item2_num = as.numeric(WFPB_micronutrientsOfConcern_folate) - 1,
         item3_num = as.numeric(WFPB_micronutrientsOfConcern_choline) - 1,
         item4_num = as.numeric(WFPB_micronutrientsOfConcern_thiamine) - 1,
         item5_num = as.numeric(WFPB_micronutrientsOfConcern_potassium) - 1,
         item6_num = as.numeric(WFPB_micronutrientsOfConcern_SComega3) - 1) %>% 
  rowwise() %>% 
  mutate(knowledge_score_uncriticalNutrients = 1 - mean(c(item1_num, item2_num, item3_num, item4_num, item5_num, item6_num))) %>% 
  ungroup()

# boxplot
ggplot(dat, aes(x = "", y = knowledge_score_uncriticalNutrients)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Uncritical micronutrients knowledge score") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Code
# alternative barplot
ggplot(dat, aes(x = knowledge_score_uncriticalNutrients, y = after_stat(count / sum(count)))) +
  geom_bar(fill = col_vector["knowledge"]) +
  scale_y_continuous("relative frequency", labels = scales::label_percent()) +
  scale_x_continuous("uncritical micronutrients knowledge score",
                     breaks = (0:6)/6, labels = paste0(0:6, "/6"),
                     limits = c(0, NA)) +
  theme(panel.grid.minor = element_blank())

Univariate description

Some general statistics:

  • arithmetic mean score (min – max): 78% (33% – 100%)
  • quantiles: 50% (2.5% quantile), 50% (5% quantile), 67% (25% quantile), 83% (50% quantile), 83% (75% quantile), 100% (95% quantile), 100% (97.5% quantile)
  • 14% of the respondents scored 100%
  • 64% of the respondents scored 75% or higher
  • 100% of the respondents scored 50% or higher
  • 100% of the respondents scored 25% or higher
  • 0% of the respondents scored 0%

The following plots are conditioned on subspecialty, age, years of practice

Code
plotDat_list <- lapply(1:length(subspecialties_varVector), function(i) {

  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # rename the respective subspecialty column, for easier handling
  dat_subsp <- dat
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  plot_dat <- data.frame(subspecialty              = subsp_label,
                         knowledge_score_uncriticalNutrients = dat_subsp$knowledge_score_uncriticalNutrients)
  
  return(plot_dat)
})
plot_dat <- plotDat_list %>% bind_rows() %>% 
  mutate(subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)))

ggplot(plot_dat, aes(x = subspecialty, y = knowledge_score_uncriticalNutrients)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Uncritical micronutrients knowledge score\nby subspecialty") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Code
dat %>% 
  filter(!is.na(age_cat)) %>% 
  ggplot(aes(x = age_cat, y = knowledge_score_uncriticalNutrients)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Uncritical micronutrients knowledge score\nby age") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Code
dat %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = knowledge_score_uncriticalNutrients)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Uncritical micronutrients knowledge score\nby years of practice") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Model-based change

Estimation of a Gaussian regression model:

  • dependent variable: Uncritical micronutrients knowledge score (quasi-metric, assumed as metric here)
  • independent variables: subspecialty (one dummy variable for each of the six subspecialties under focus), years of practice (in 10-year-units), education quality on WFPB diets (as numeric score with 0 = ‘No training’, 1 = ‘Very poor’ and 5 = ‘Excellent’)
  • further control variables: education
  • unused control variables: age (since strongly correlated with years of practice), gender (since dietitians are nearly exclusively female), personal diet (since the causality between a dietitian’s personal diet and the response variable is unclear)
Code
# model estimation
model <- gam(knowledge_score_uncriticalNutrients ~ 
               dieteticsSubspecialty_obesityOrWeightManagement +
               dieteticsSubspecialty_gastroenterology +
               dieteticsSubspecialty_diabetes +
               dieteticsSubspecialty_paediatrics +
               dieteticsSubspecialty_elderly +
               dieteticsSubspecialty_oncology +
               dieteticsSubspecialty_psychEating +
               dietitian_forHowLong_10yearUnits +
               education +
               education_WFPBeducation_num,
             data = dat)

# prepare the results table
results_tab_UNscore <- summary(model)$p.table %>% 
  as.data.frame() %>% 
  mutate(parameter = row.names(.)) %>% 
  rename(estimate = "Estimate",
         se       = "Std. Error",
         pvalue   = "Pr(>|t|)") %>% 
  mutate(CI_lower = estimate - qnorm(.975) * se,
         CI_upper = estimate + qnorm(.975) * se,
         pvalue   = case_when(pvalue < .0001 ~ "<.0001",
                              TRUE           ~ as.character(round(pvalue, 4))),
         param_group = case_when(grepl("Subspecialty", parameter)  ~ "area of specialty",
                                 grepl("forHowLong", parameter)    ~ "years of practice",
                                 grepl("WFPBeducation", parameter) ~ "WFPB education",
                                 grepl("education", parameter)     ~ "education"),
         param_group = factor(param_group, levels = c("area of specialty", "years of practice", "education", "WFPB education")),
         parameter    = gsub(parameter, pattern = "dieteticsSubspecialty_", replacement = "")) %>% 
  select(param_group, parameter, estimate, CI_lower, CI_upper, pvalue)
row.names(results_tab_UNscore) <- NULL
results_tab_UNscore %>% 
  mutate(estimate = round(estimate, 4),
         CI_lower = round(CI_lower, 4),
         CI_upper = round(CI_upper, 4)) %>% 
  kable() %>% 
  kable_styling()
param_group parameter estimate CI_lower CI_upper pvalue
NA (Intercept) 0.7410 0.6988 0.7832 <.0001
area of specialty obesityOrWeightManagementyes 0.0012 -0.0355 0.0380 0.948
area of specialty gastroenterologyyes 0.0015 -0.0383 0.0414 0.9396
area of specialty diabetesyes -0.0090 -0.0515 0.0335 0.6788
area of specialty paediatricsyes 0.0145 -0.0260 0.0550 0.4845
area of specialty elderlyyes 0.0062 -0.0361 0.0485 0.7742
area of specialty oncologyyes 0.0016 -0.0469 0.0501 0.9483
area of specialty psychEatingyes 0.0086 -0.0404 0.0576 0.7305
years of practice dietitian_forHowLong_10yearUnits 0.0160 0.0005 0.0314 0.0441
education educationPostgraduate degree 0.0324 0.0015 0.0633 0.0409
education educationPhD 0.0393 -0.0325 0.1111 0.2845
WFPB education education_WFPBeducation_num -0.0002 -0.0143 0.0138 0.9738
Code
# effect plot
results_tab_UNscore %>% 
  filter(parameter != "(Intercept)") %>% 
  ggplot(aes(x = parameter)) +
  geom_hline(yintercept = 0, lty = 2, col = "gray50") +
  geom_pointrange(aes(y = estimate, ymin = CI_lower, ymax = CI_upper, col = param_group)) +
  facet_wrap(~ param_group, scales = "free_x", nrow = 1) +
  theme(axis.text.x     = element_text(angle = 45, hjust = 1),
        axis.title.x    = element_blank(),
        legend.position = "none")

R²: -0.001
Share of explained deviance: 3.2%

Residual structure looks acceptable.

Code
# # check the residuals
# resid_dat <- data.frame(prediction = model$fitted.values,
#                         residuum   = model$residuals)
# ggplot(resid_dat, aes(x = prediction, y = residuum)) +
#   geom_point() +
#   geom_smooth(se = FALSE) +
#   ggtitle("Model residuals vs. predicted values")

WFPB: Improved diseases knowledge score

Does knowledge about which disease’s conditions improve through WFPB diets change or is determined by i) age or ii) years practising as a dietitian?

We define a improved diseases WFPB knowledge score based on these 13 diseases:

  1. Heart disease
  2. High cholesterol
  3. Hypertension
  4. Diabetes
  5. Obesity
  6. Stroke
  7. Chronic kidney disease
  8. Cancers
  9. Fatty liver disease
  10. Alzeimers dementia
  11. Vascular Dementia
  12. Depression
  13. Inflammatory disease

Based on these 13 diseases, the knowledge score for every person is calculated by taking the average of a person’s responses to these items. The obtained score then lies between zero and one and reflects the share of the above diseases which were correctly answered as benefiting by a WFPB diet through improved conditions.

Code
dat <- dat %>% 
  mutate(item1_num = as.numeric(Pbdiet_improvedConditions_heartDisease) - 1,
         item2_num = as.numeric(Pbdiet_improvedConditions_highCholesterol) - 1,
         item3_num = as.numeric(Pbdiet_improvedConditions_hypertension) - 1,
         item4_num = as.numeric(Pbdiet_improvedConditions_T2Diabetes) - 1,
         item5_num = as.numeric(Pbdiet_improvedConditions_obesity) - 1,
         item6_num = as.numeric(Pbdiet_improvedConditions_stroke) - 1,
         item7_num = as.numeric(Pbdiet_improvedConditions_chronicKidney) - 1,
         item8_num = as.numeric(Pbdiet_improvedConditions_certainCancers) - 1,
         item9_num = as.numeric(Pbdiet_improvedConditions_fattyLiver) - 1,
         item10_num = as.numeric(Pbdiet_improvedConditions_alzheimer) - 1,
         item11_num = as.numeric(Pbdiet_improvedConditions_vascularDementia) - 1,
         item12_num = as.numeric(Pbdiet_improvedConditions_depression) - 1,
         item13_num = as.numeric(Pbdiet_improvedConditions_IBD) - 1) %>% 
  rowwise() %>% 
  mutate(knowledge_score_diseases = mean(c(item1_num, item2_num, item3_num, item4_num, item5_num, item6_num, item7_num, item8_num, item9_num, item10_num, item11_num, item12_num, item13_num))) %>% 
  ungroup()

# boxplot
ggplot(dat, aes(x = "", y = knowledge_score_diseases)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Improved diseases knowledge score") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Code
# alternative barplot
ggplot(dat, aes(x = knowledge_score_diseases, y = after_stat(count / sum(count)))) +
  geom_bar(fill = col_vector["knowledge"]) +
  scale_y_continuous("relative frequency", labels = scales::label_percent()) +
  scale_x_continuous("improved diseases knowledge score", labels = scales::label_percent()) +
  theme(panel.grid.minor = element_blank())

Univariate description

Some general statistics:

  • arithmetic mean score (min – max): 54% (0% – 100%)
  • quantiles: 9% (2.5% quantile), 15% (5% quantile), 38% (25% quantile), 54% (50% quantile), 69% (75% quantile), 100% (95% quantile), 100% (97.5% quantile)
  • 7% of the respondents scored 100%
  • 22% of the respondents scored 75% or higher
  • 55% of the respondents scored 50% or higher
  • 84% of the respondents scored 25% or higher
  • 1% of the respondents scored 0%

The following plots are conditioned on subspecialty, age, years of practice

Code
plotDat_list <- lapply(1:length(subspecialties_varVector), function(i) {

  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # rename the respective subspecialty column, for easier handling
  dat_subsp <- dat
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  plot_dat <- data.frame(subspecialty             = subsp_label,
                         knowledge_score_diseases = dat_subsp$knowledge_score_diseases)
  
  return(plot_dat)
})
plot_dat <- plotDat_list %>% bind_rows() %>% 
  mutate(subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)))

ggplot(plot_dat, aes(x = subspecialty, y = knowledge_score_diseases)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Improved diseases knowledge score\nby subspecialty") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Code
dat %>% 
  filter(!is.na(age_cat)) %>% 
  ggplot(aes(x = age_cat, y = knowledge_score_diseases)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Improved diseases knowledge score\nby age") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Code
dat %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = knowledge_score_diseases)) +
  geom_boxplot(color = col_vector["knowledge"]) +
  ggtitle("WFPB: Improved diseases knowledge score\nby years of practice") +
  ylab("knowledge score") +
  theme(axis.title.x       = element_blank(),
        panel.grid.major.x = element_blank())

Model-based change

Estimation of a Gaussian regression model:

  • dependent variable: Diseases knowledge score (quasi-metric, assumed as metric here)
  • independent variables: subspecialty (one dummy variable for each of the six subspecialties under focus), years of practice (in 10-year-units), education quality on WFPB diets (as numeric score with 0 = ‘No training’, 1 = ‘Very poor’ and 5 = ‘Excellent’)
  • further control variables: education
  • unused control variables: age (since strongly correlated with years of practice), gender (since dietitians are nearly exclusively female), personal diet (since the causality between a dietitian’s personal diet and the response variable is unclear)
Code
# model estimation
model <- gam(knowledge_score_diseases ~ 
               dieteticsSubspecialty_obesityOrWeightManagement +
               dieteticsSubspecialty_gastroenterology +
               dieteticsSubspecialty_diabetes +
               dieteticsSubspecialty_paediatrics +
               dieteticsSubspecialty_elderly +
               dieteticsSubspecialty_oncology +
               dieteticsSubspecialty_psychEating +
               dietitian_forHowLong_10yearUnits +
               education +
               education_WFPBeducation_num,
             data = dat)

# prepare the results table
results_tab_IDscore <- summary(model)$p.table %>% 
  as.data.frame() %>% 
  mutate(parameter = row.names(.)) %>% 
  rename(estimate = "Estimate",
         se       = "Std. Error",
         pvalue   = "Pr(>|t|)") %>% 
  mutate(CI_lower = estimate - qnorm(.975) * se,
         CI_upper = estimate + qnorm(.975) * se,
         pvalue   = case_when(pvalue < .0001 ~ "<.0001",
                              TRUE           ~ as.character(round(pvalue, 4))),
         param_group = case_when(grepl("Subspecialty", parameter)  ~ "area of specialty",
                                 grepl("forHowLong", parameter)    ~ "years of practice",
                                 grepl("WFPBeducation", parameter) ~ "WFPB education",
                                 grepl("education", parameter)     ~ "education"),
         param_group = factor(param_group, levels = c("area of specialty", "years of practice", "education", "WFPB education")),
         parameter    = gsub(parameter, pattern = "dieteticsSubspecialty_", replacement = "")) %>% 
  select(param_group, parameter, estimate, CI_lower, CI_upper, pvalue)
row.names(results_tab_IDscore) <- NULL
results_tab_IDscore %>% 
  mutate(estimate = round(estimate, 4),
         CI_lower = round(CI_lower, 4),
         CI_upper = round(CI_upper, 4)) %>% 
  kable() %>% 
  kable_styling()
param_group parameter estimate CI_lower CI_upper pvalue
NA (Intercept) 0.5559 0.4808 0.6310 <.0001
area of specialty obesityOrWeightManagementyes 0.0368 -0.0285 0.1022 0.2699
area of specialty gastroenterologyyes 0.0421 -0.0287 0.1129 0.2448
area of specialty diabetesyes 0.0248 -0.0507 0.1004 0.52
area of specialty paediatricsyes 0.0647 -0.0073 0.1368 0.0791
area of specialty elderlyyes -0.0083 -0.0835 0.0669 0.8286
area of specialty oncologyyes -0.0400 -0.1262 0.0461 0.3632
area of specialty psychEatingyes 0.0685 -0.0186 0.1556 0.1243
years of practice dietitian_forHowLong_10yearUnits -0.0163 -0.0438 0.0112 0.2456
education educationPostgraduate degree 0.0441 -0.0110 0.0991 0.1175
education educationPhD -0.0028 -0.1305 0.1249 0.9655
WFPB education education_WFPBeducation_num -0.0308 -0.0558 -0.0058 0.0165
Code
# effect plot
results_tab_IDscore %>% 
  filter(parameter != "(Intercept)") %>% 
  ggplot(aes(x = parameter)) +
  geom_hline(yintercept = 0, lty = 2, col = "gray50") +
  geom_pointrange(aes(y = estimate, ymin = CI_lower, ymax = CI_upper, col = param_group)) +
  facet_wrap(~ param_group, scales = "free_x", nrow = 1) +
  theme(axis.text.x     = element_text(angle = 45, hjust = 1),
        axis.title.x    = element_blank(),
        legend.position = "none")

R²: 0.023
Share of explained deviance: 5.6%

Residual structure looks acceptable.

Code
# # check the residuals
# resid_dat <- data.frame(prediction = model$fitted.values,
#                         residuum   = model$residuals)
# ggplot(resid_dat, aes(x = prediction, y = residuum)) +
#   geom_point() +
#   geom_smooth(se = FALSE) +
#   ggtitle("Model residuals vs. predicted values")

Beliefs research questions

WFPB: Long-term sustainability

Do they believe WFPB are sustainable long-term?

Code
dat %>% 
  ggplot(aes(y = "", fill = WFPB_isSustainableLongTerm)) +
  geom_bar(position = "fill") +
  scale_x_continuous("Relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "BrBG") +
  ggtitle("Long-term sustainability of WFPB diets") +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

Univariate description

by subspecialty, age, years of practice, and personal dietary pattern.

Code
# prepare the dataset
plot_dat_list <- lapply(1:length(subspecialties_varVector), function(i) {
  
  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # filter the dataset for dietitians with the specific subspecialty
  dat_subsp <- dat
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  results_dat <- data.frame(WFPB_isSustainableLongTerm = levels(dat_subsp$WFPB_isSustainableLongTerm),
             share_yes                   = prop.table(table(dat_subsp$WFPB_isSustainableLongTerm)) %>% as.vector,
             subspecialty                = subsp_label) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat_sustainability_sub <- plot_dat_list %>% bind_rows() %>% 
  mutate(subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)),
         WFPB_isSustainableLongTerm = factor(WFPB_isSustainableLongTerm, levels = levels(dat$WFPB_isSustainableLongTerm)))

# heatmap
plot_dat_sustainability_sub %>% 
  ggplot(aes(x = subspecialty, y = WFPB_isSustainableLongTerm, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["beliefs"]),
                      labels = scales::label_percent()) +
  ggtitle("Long-term sustainability of WFPB diets\nby subspecialties") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Code
# prepare the dataset
plot_dat_list <- lapply(levels(dat$age_cat), function(x) {
  
  # filter the dataset for dietitians with the specific age group
  dat_age <- dat %>% filter(age_cat == x)
  
  results_dat <- data.frame(WFPB_isSustainableLongTerm = levels(dat_age$WFPB_isSustainableLongTerm),
                            share_yes                   = prop.table(table(dat_age$WFPB_isSustainableLongTerm)) %>% as.vector,
                            age_cat                     = x) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat_sustainability_age <- plot_dat_list %>% bind_rows() %>% 
  mutate(age_cat                     = factor(age_cat, levels = levels(dat$age_cat)),
         WFPB_isSustainableLongTerm  = factor(WFPB_isSustainableLongTerm, levels = levels(dat$WFPB_isSustainableLongTerm)))

# heatmap
plot_dat_sustainability_age %>% 
  ggplot(aes(x = age_cat, y = WFPB_isSustainableLongTerm, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["beliefs"]),
                      labels = scales::label_percent()) +
  ggtitle("Long-term sustainability of WFPB diets\nby age groups") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Code
# prepare the dataset
plot_dat_list <- lapply(levels(dat$dietitian_forHowLong_cat), function(x) {
  
  # filter the dataset for dietitians with the specific years of practice
  dat_practice <- dat %>% filter(dietitian_forHowLong_cat == x)
  
  results_dat <- data.frame(WFPB_isSustainableLongTerm  = levels(dat_practice$WFPB_isSustainableLongTerm),
                            share_yes                   = prop.table(table(dat_practice$WFPB_isSustainableLongTerm)) %>% as.vector,
                            dietitian_forHowLong_cat    = x) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat_sustainability_years <- plot_dat_list %>% bind_rows() %>% 
  mutate(dietitian_forHowLong_cat    = factor(dietitian_forHowLong_cat, levels = levels(dat$dietitian_forHowLong_cat)),
         WFPB_isSustainableLongTerm  = factor(WFPB_isSustainableLongTerm, levels = levels(dat$WFPB_isSustainableLongTerm)))

# heatmap
plot_dat_sustainability_years %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = WFPB_isSustainableLongTerm, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["beliefs"]),
                      labels = scales::label_percent()) +
  ggtitle("Long-term sustainability of WFPB diets\nby dietitians' years of practice") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Code
# prepare the dataset
plot_dat_list <- lapply(levels(dat$diet_personal), function(x) {
  
  # filter the dataset for dietitians with the specific diet
  dat_diet <- dat %>% filter(diet_personal == x)
  
  results_dat <- data.frame(WFPB_isSustainableLongTerm  = levels(dat_diet$WFPB_isSustainableLongTerm),
                            share_yes                   = prop.table(table(dat_diet$WFPB_isSustainableLongTerm)) %>% as.vector,
                            diet_personal               = x) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat <- plot_dat_list %>% bind_rows() %>% 
  mutate(diet_personal              = factor(diet_personal, levels = levels(dat$diet_personal)),
         WFPB_isSustainableLongTerm = factor(WFPB_isSustainableLongTerm, levels = levels(dat$WFPB_isSustainableLongTerm)))

# heatmap
plot_dat %>% 
  ggplot(aes(x = diet_personal, y = WFPB_isSustainableLongTerm, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["beliefs"]),
                      labels = scales::label_percent()) +
  ggtitle("Long-term sustainability of WFPB diets\nby dietitians' personal diet") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Model-based change

Estimation of a Logistic regression model:

  • dependent variable: binarized long-term sustainability of WFPB diets (‘(Strongly) Agree’ vs. ‘Not sure / (Strongly) Disagree’)
  • independent variables: subspecialty (one dummy variable for each of the six subspecialties under focus), years of practice (in 10-year-units), education quality on WFPB diets (as numeric score with 0 = ‘No training’, 1 = ‘Very poor’ and 5 = ‘Excellent’)
  • further control variables: education
  • unused control variables: age (since strongly correlated with years of practice), gender (since dietitians are nearly exclusively female), personal diet (since the causality between a dietitian’s personal diet and the response variable is unclear)
Code
# calculate the response variable
dat <- dat %>% 
  mutate(WFPBsustainability_binary = case_when(WFPB_isSustainableLongTerm %in% c("Agree","Strongly agree")                  ~ "(Strongly) Agree",
                                               WFPB_isSustainableLongTerm %in% c("Not sure","Disagree","Strongly disagree") ~ "Not sure / (Strongly) Disagree",
                                               TRUE                                                                         ~ NA_character_),
         WFPBsustainability_binary = factor(WFPBsustainability_binary, levels = c("Not sure / (Strongly) Disagree","(Strongly) Agree")))

# model estimation
model <- gam(WFPBsustainability_binary ~ 
               dieteticsSubspecialty_obesityOrWeightManagement +
               dieteticsSubspecialty_gastroenterology +
               dieteticsSubspecialty_diabetes +
               dieteticsSubspecialty_paediatrics +
               dieteticsSubspecialty_elderly +
               dieteticsSubspecialty_oncology +
               dieteticsSubspecialty_psychEating +
               dietitian_forHowLong_10yearUnits +
               education +
               education_WFPBeducation_num,
             family = binomial(link = "logit"),
             data = dat)

# prepare the results table
results_tab_longTermSus <- summary(model)$p.table %>% 
  as.data.frame() %>% 
  mutate(parameter = row.names(.)) %>% 
  rename(estimate = "Estimate",
         se       = "Std. Error",
         pvalue   = "Pr(>|z|)") %>% 
  mutate(CI_lower     = estimate - qnorm(.975) * se,
         CI_upper     = estimate + qnorm(.975) * se,
         exp_estimate = exp(estimate),
         exp_CIlower  = exp(CI_lower),
         exp_CIupper  = exp(CI_upper),
         pvalue       = case_when(pvalue < .0001 ~ "<.0001",
                                  TRUE           ~ as.character(round(pvalue, 4))),
         param_group  = case_when(grepl("Subspecialty", parameter)  ~ "area of specialty",
                                  grepl("forHowLong", parameter)    ~ "years of practice",
                                  grepl("WFPBeducation", parameter) ~ "WFPB education",
                                  grepl("education", parameter)     ~ "education"),
         param_group  = factor(param_group, levels = c("area of specialty", "years of practice", "education", "WFPB education")),
         parameter    = gsub(parameter, pattern = "dieteticsSubspecialty_", replacement = "")) %>% 
  select(param_group, parameter, exp_estimate, exp_CIlower, exp_CIupper, pvalue)
row.names(results_tab_longTermSus) <- NULL
results_tab_longTermSus %>% 
  mutate(exp_estimate = round(exp_estimate, 4),
         exp_CIlower  = round(exp_CIlower, 4),
         exp_CIupper  = round(exp_CIupper, 4)) %>% 
  kable() %>% 
  kable_styling()
param_group parameter exp_estimate exp_CIlower exp_CIupper pvalue
NA (Intercept) 1.6900 0.8630 3.3095 0.126
area of specialty obesityOrWeightManagementyes 1.0335 0.5723 1.8665 0.913
area of specialty gastroenterologyyes 1.6959 0.8613 3.3396 0.1265
area of specialty diabetesyes 1.1104 0.5601 2.2014 0.7641
area of specialty paediatricsyes 1.2391 0.6390 2.4027 0.5258
area of specialty elderlyyes 1.2657 0.6332 2.5304 0.5049
area of specialty oncologyyes 0.3920 0.1891 0.8127 0.0118
area of specialty psychEatingyes 1.0285 0.4706 2.2477 0.9438
years of practice dietitian_forHowLong_10yearUnits 0.9682 0.7564 1.2392 0.7973
education educationPostgraduate degree 1.7272 1.0538 2.8310 0.0302
education educationPhD 1.5183 0.4848 4.7546 0.4734
WFPB education education_WFPBeducation_num 0.9328 0.7436 1.1700 0.5472
Code
# effect plot
results_tab_longTermSus %>% 
  filter(parameter != "(Intercept)") %>% 
  ggplot(aes(x = parameter)) +
  geom_hline(yintercept = 1, lty = 2, col = "gray50") +
  geom_pointrange(aes(y = exp_estimate, ymin = exp_CIlower, ymax = exp_CIupper, col = param_group)) +
  scale_y_continuous("Odds Ratio, on log2-transformed scale", trans = "log2") +
  facet_wrap(~ param_group, scales = "free_x", nrow = 1) +
  theme(axis.text.x     = element_text(angle = 45, hjust = 1),
        axis.title.x    = element_blank(),
        legend.position = "none")

The model’s goodness-of-fit is determined by calculating the AUC (Area Under the Curve). Therefore, the dataset is first randomly split into a training dataset (80% of the data) and a test dataset (20% of the data). The regression model is then re-estimated on the training data, and the AUC is calculated on the test data.

Code
# basis: split data into a training set (80% of the data) and a test set (20%)
set.seed(2024)
train_obs <- sample(1:nrow(dat), size = round(0.8*nrow(dat)), replace = FALSE)
dat_train <- dat %>% slice(train_obs)
dat_test  <- dat %>% slice(-train_obs)
Code
# 1) estimate the model on the training set
model_train <- gam(WFPBsustainability_binary ~ 
                     dieteticsSubspecialty_obesityOrWeightManagement +
                     dieteticsSubspecialty_gastroenterology +
                     dieteticsSubspecialty_diabetes +
                     dieteticsSubspecialty_paediatrics +
                     dieteticsSubspecialty_elderly +
                     dieteticsSubspecialty_oncology +
                     dieteticsSubspecialty_psychEating +
                     dietitian_forHowLong_10yearUnits +
                     education,
                   family = binomial(link = "logit"),
                   data = dat_train)
# 2) calculate ROC curve and the corresponding AUC value
pred <- data.frame("outcome" = as.integer(dat_test$WFPBsustainability_binary) - 1,
                   "pred"    = predict(model_train, newdata = dat_test))
gg <- ggplot(pred, aes(m = pred, d = outcome)) + geom_roc() + ggtitle("ROC curve")
auc <- calc_auc(gg)$AUC

AUC: 0.55

WFPB: Long-term adherence under co-morbidities

Do they believe clients with co-morbidities would adhere long-term to a WFPB diet?

Code
dat %>% 
  ggplot(aes(y = "", fill = ChronicConditions_longTermWFPBadherence)) +
  geom_bar(position = "fill") +
  scale_x_continuous("Relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "BrBG") +
  ggtitle("Long-term adherence under co-morbidities") +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

WFPB: Clients’ main motivation

What is the main motivation they have come across with their clients in adopting a WFPB diet?

Code
dat %>% 
  select(starts_with("Pbdiet_primaryMotivationClients")) %>% 
  select(-Pbdiet_primaryMotivationClients) %>% 
  pivot_longer(cols = everything()) %>% 
  group_by(name) %>% 
  mutate(sum_yes = sum(value == "yes")) %>% 
  ungroup() %>% 
  arrange(desc(sum_yes)) %>% 
  mutate(name = gsub(pattern = "Pbdiet_primaryMotivationClients", replacement = "", name),
         name = factor(name, levels = rev(unique(name)))) %>% 
  ggplot(aes(y = name, fill = value)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("lightgray", unname(col_vector["beliefs"]))) +
  scale_x_continuous(labels = scales::label_percent()) +
  xlab("Agreement [%]") +
  ggtitle("Clients' primary motivation for plant-based diet") +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

Regularity of WFPB recommendation

How often would they recommend a WFPB diet in clinical practice?

Excluding 13 people from the following plots with value Not applicable.

Code
dat %>% 
  filter(WFPBrecommendation_howOften != "Not applicable") %>% 
  group_by(WFPBrecommendation_howOften) %>% 
  summarize(n = n()) %>% 
  mutate(rel_freq = n / sum(n)) %>% 
  ggplot(aes(x = WFPBrecommendation_howOften, weight = rel_freq)) +
  geom_bar(fill = col_vector["beliefs"]) +
  scale_y_continuous(labels = scales::label_percent()) +
  ggtitle("How often would you recommend a WFPB diet to your patients?") +
  theme(axis.title = element_blank())

Univariate description

by subspecialty, age, years of practice.

Code
# prepare the dataset
plot_dat_list <- lapply(1:length(subspecialties_varVector), function(i) {
  
  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # filter the dataset for dietitians with the specific subspecialty
  dat_subsp <- dat %>% filter(WFPBrecommendation_howOften != "Not applicable") %>% 
    mutate(WFPBrecommendation_howOften = droplevels(WFPBrecommendation_howOften))
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  results_dat <- data.frame(WFPBrecommendation_howOften = levels(dat_subsp$WFPBrecommendation_howOften),
             share_yes                   = prop.table(table(dat_subsp$WFPBrecommendation_howOften)) %>% as.vector,
             subspecialty                = subsp_label) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat <- plot_dat_list %>% bind_rows() %>% 
  mutate(subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)),
         WFPBrecommendation_howOften = factor(WFPBrecommendation_howOften, levels = levels(dat$WFPBrecommendation_howOften)))

# heatmap
plot_dat %>% 
  ggplot(aes(x = subspecialty, y = WFPBrecommendation_howOften, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["beliefs"]),
                      labels = scales::label_percent()) +
  ggtitle("Regularity of WFPB recommendation\nby subspecialties") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Code
# prepare the dataset
plot_dat_list <- lapply(levels(dat$age_cat), function(x) {
  
  # filter the dataset for dietitians with the specific age group
  dat_age <- dat %>% 
    filter(WFPBrecommendation_howOften != "Not applicable") %>% 
    mutate(WFPBrecommendation_howOften = droplevels(WFPBrecommendation_howOften)) %>% 
    filter(age_cat == x)
  
  results_dat <- data.frame(WFPBrecommendation_howOften = levels(dat_age$WFPBrecommendation_howOften),
                            share_yes                   = prop.table(table(dat_age$WFPBrecommendation_howOften)) %>% as.vector,
                            age_cat                     = x) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat <- plot_dat_list %>% bind_rows() %>% 
  mutate(age_cat                     = factor(age_cat, levels = levels(dat$age_cat)),
         WFPBrecommendation_howOften = factor(WFPBrecommendation_howOften, levels = levels(dat$WFPBrecommendation_howOften)))

# heatmap
plot_dat %>% 
  ggplot(aes(x = age_cat, y = WFPBrecommendation_howOften, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["beliefs"]),
                      labels = scales::label_percent()) +
  ggtitle("Regularity of WFPB recommendation\nby age groups") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Code
# prepare the dataset
plot_dat_list <- lapply(levels(dat$dietitian_forHowLong_cat), function(x) {
  
  # filter the dataset for dietitians with the specific years of practice
  dat_practice <- dat %>% 
    filter(WFPBrecommendation_howOften != "Not applicable") %>% 
    mutate(WFPBrecommendation_howOften = droplevels(WFPBrecommendation_howOften)) %>% 
    filter(dietitian_forHowLong_cat == x)
  
  results_dat <- data.frame(WFPBrecommendation_howOften = levels(dat_practice$WFPBrecommendation_howOften),
                            share_yes                   = prop.table(table(dat_practice$WFPBrecommendation_howOften)) %>% as.vector,
                            dietitian_forHowLong_cat    = x) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat <- plot_dat_list %>% bind_rows() %>% 
  mutate(dietitian_forHowLong_cat    = factor(dietitian_forHowLong_cat, levels = levels(dat$dietitian_forHowLong_cat)),
         WFPBrecommendation_howOften = factor(WFPBrecommendation_howOften, levels = levels(dat$WFPBrecommendation_howOften)))

# heatmap
plot_dat %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = WFPBrecommendation_howOften, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["beliefs"]),
                      labels = scales::label_percent()) +
  ggtitle("Regularity of WFPB recommendation\nby dietitians' years of practice") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Model-based change

Estimation of a Logistic regression model:

  • dependent variable: binarized regularity of WFPB recommendation (‘Sometimes / Often / Always’ vs. ‘Rarely / Never’)
  • independent variables: subspecialty (one dummy variable for each of the six subspecialties under focus), years of practice (in 10-year-units), education quality on WFPB diets (as numeric score with 0 = ‘No training’, 1 = ‘Very poor’ and 5 = ‘Excellent’)
  • further control variables: education, educational resources on WFPB diets (as numeric score with 0 = ‘Strongly disagree’ and 4 = ‘Strongly agree’), workplace support on advocating for WFPB diets (as numeric score with 0 = ‘Strongly unsupported’ and 4 = ‘Strongly supported’)
  • unused control variables: age (since strongly correlated with years of practice), gender (since dietitians are nearly exclusively female), personal diet (since the causality between a dietitian’s personal diet and the response variable is unclear)
Code
# calculate the response variable
dat <- dat %>% 
  mutate(WFPBrecommendation_binary = case_when(WFPBrecommendation_howOften %in% c("Sometimes","Often","Always") ~ "Sometimes / Often / Always",
                                               WFPBrecommendation_howOften %in% c("Rarely","Never")             ~ "Rarely / Never",
                                               TRUE                                                             ~ NA_character_),
         WFPBrecommendation_binary = factor(WFPBrecommendation_binary, levels = c("Rarely / Never", "Sometimes / Often / Always")))

# make numeric scores out of some control variables
dat <- dat %>% 
  mutate(PBdiet_enoughEduResources_num        = as.numeric(Pbdiet_enoughProperEducationalResources) - 1,
         workplace_supportForWFPBadvocacy_num = as.numeric(workplace_supportForWFPBadvocacy) - 1)

# model estimation
model <- gam(WFPBrecommendation_binary ~ 
               dieteticsSubspecialty_obesityOrWeightManagement +
               dieteticsSubspecialty_gastroenterology +
               dieteticsSubspecialty_diabetes +
               dieteticsSubspecialty_paediatrics +
               dieteticsSubspecialty_elderly +
               dieteticsSubspecialty_oncology +
               dieteticsSubspecialty_psychEating +
               dietitian_forHowLong_10yearUnits +
               education +
               education_WFPBeducation_num +
               PBdiet_enoughEduResources_num +
               workplace_supportForWFPBadvocacy_num,
             family = binomial(link = "logit"),
             data = dat)

# prepare the results table
results_tab <- summary(model)$p.table %>% 
  as.data.frame() %>% 
  mutate(parameter = row.names(.)) %>% 
  rename(estimate = "Estimate",
         se       = "Std. Error",
         pvalue   = "Pr(>|z|)") %>% 
  mutate(CI_lower     = estimate - qnorm(.975) * se,
         CI_upper     = estimate + qnorm(.975) * se,
         exp_estimate = exp(estimate),
         exp_CIlower  = exp(CI_lower),
         exp_CIupper  = exp(CI_upper),
         pvalue       = case_when(pvalue < .0001 ~ "<.0001",
                                  TRUE           ~ as.character(round(pvalue, 4))),
         param_group  = case_when(grepl("Subspecialty", parameter) ~ "area of specialty",
                                  grepl("forHowLong", parameter)   ~ "years of practice",
                                  grepl("_num", parameter)         ~ "barriers",
                                  grepl("education", parameter)    ~ "education"),
         param_group  = factor(param_group, levels = c("
                                                       ", "years of practice", "education", "barriers")),
         parameter    = gsub(parameter, pattern = "dieteticsSubspecialty_", replacement = "")) %>% 
  select(param_group, parameter, exp_estimate, exp_CIlower, exp_CIupper, pvalue)
row.names(results_tab) <- NULL
results_tab %>% 
  mutate(exp_estimate = round(exp_estimate, 4),
         exp_CIlower  = round(exp_CIlower, 4),
         exp_CIupper  = round(exp_CIupper, 4)) %>% 
  kable() %>% 
  kable_styling()
param_group parameter exp_estimate exp_CIlower exp_CIupper pvalue
NA (Intercept) 0.1167 0.0469 0.2901 <.0001
NA obesityOrWeightManagementyes 1.8364 1.0206 3.3044 0.0426
NA gastroenterologyyes 1.3231 0.6915 2.5314 0.3977
NA diabetesyes 1.1655 0.5967 2.2765 0.6539
NA paediatricsyes 0.7996 0.4083 1.5657 0.5142
NA elderlyyes 0.9868 0.4907 1.9845 0.9702
NA oncologyyes 0.5132 0.2211 1.1913 0.1205
NA psychEatingyes 1.2272 0.5556 2.7105 0.6126
years of practice dietitian_forHowLong_10yearUnits 1.2509 0.9771 1.6013 0.0757
education educationPostgraduate degree 2.4540 1.4715 4.0926 6e-04
education educationPhD 4.4036 1.4067 13.7854 0.0109
barriers education_WFPBeducation_num 0.9333 0.7401 1.1768 0.5594
barriers PBdiet_enoughEduResources_num 1.1198 0.8532 1.4697 0.4147
barriers workplace_supportForWFPBadvocacy_num 1.3925 1.0849 1.7874 0.0093
Code
# effect plot
results_tab %>% 
  filter(parameter != "(Intercept)") %>% 
  filter(exp_CIupper < 1000) %>% # hide the 'WFPBeducation = excellent' subgroup as of only one observation in it
  ggplot(aes(x = parameter)) +
  geom_hline(yintercept = 1, lty = 2, col = "gray50") +
  geom_pointrange(aes(y = exp_estimate, ymin = exp_CIlower, ymax = exp_CIupper, col = param_group)) +
  scale_y_continuous("Odds Ratio, on log2-transformed scale", trans = "log2") +
  facet_wrap(~ param_group, scales = "free_x", nrow = 1) +
  theme(axis.text.x     = element_text(angle = 45, hjust = 1),
        axis.title.x    = element_blank(),
        legend.position = "none")

The model’s goodness-of-fit is determined by calculating the AUC (Area Under the Curve). Therefore, the dataset is first randomly split into a training dataset (80% of the data) and a test dataset (20% of the data). The regression model is then re-estimated on the training data, and the AUC is calculated on the test data.

Code
# basis: split data into a training set (80% of the data) and a test set (20%)
set.seed(2024)
train_obs <- sample(1:nrow(dat), size = round(0.8*nrow(dat)), replace = FALSE)
dat_train <- dat %>% slice(train_obs)
dat_test  <- dat %>% slice(-train_obs)
Code
# 1) estimate the model on the training set
model_train <- gam(WFPBrecommendation_binary ~ 
                     dieteticsSubspecialty_obesityOrWeightManagement +
                     dieteticsSubspecialty_gastroenterology +
                     dieteticsSubspecialty_diabetes +
                     dieteticsSubspecialty_paediatrics +
                     dieteticsSubspecialty_elderly +
                     dieteticsSubspecialty_oncology +
                     dieteticsSubspecialty_psychEating +
                     dietitian_forHowLong_10yearUnits +
                     education,
                   family = binomial(link = "logit"),
                   data = dat_train)
# 2) calculate ROC curve and the corresponding AUC value
pred <- data.frame("outcome" = as.integer(dat_test$WFPBrecommendation_binary) - 1,
                   "pred"    = predict(model_train, newdata = dat_test))
gg <- ggplot(pred, aes(m = pred, d = outcome)) + geom_roc() + ggtitle("ROC curve")
auc <- calc_auc(gg)$AUC

AUC: 0.61

WFPB: Availability in hospitals

Do they believe that WFPB diets should be integrated into hospitals or other healthcare facilities?

Code
dat %>% 
  ggplot(aes(y = "", fill = withinHealthcare_PbdietShouldBeOption)) +
  geom_bar(position = "fill") +
  scale_x_continuous("Relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "BrBG") +
  guides(fill = guide_legend(reverse = TRUE)) +
  ggtitle("Should WFPB diets be available in hospitals?") +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

Univariate description

by subspecialty and years of practice

Code
# prepare the dataset
plot_dat_list <- lapply(1:length(subspecialties_varVector), function(i) {
  
  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # filter the dataset for dietitians with the specific subspecialty
  dat_subsp <- dat
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  results_dat <- data.frame(withinHealthcare_PbdietShouldBeOption = levels(dat_subsp$withinHealthcare_PbdietShouldBeOption),
             share_yes                   = prop.table(table(dat_subsp$withinHealthcare_PbdietShouldBeOption)) %>% as.vector,
             subspecialty                = subsp_label) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat_hospitals_sub <- plot_dat_list %>% bind_rows() %>% 
  mutate(subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)),
         withinHealthcare_PbdietShouldBeOption = factor(withinHealthcare_PbdietShouldBeOption, levels = levels(dat$withinHealthcare_PbdietShouldBeOption)))

# heatmap
plot_dat_hospitals_sub %>% 
  ggplot(aes(x = subspecialty, y = withinHealthcare_PbdietShouldBeOption, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["beliefs"]),
                      labels = scales::label_percent()) +
  ggtitle("Should WFPB diets be available in hospitals?\nby subspecialties") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Code
# prepare the dataset
plot_dat_list <- lapply(levels(dat$dietitian_forHowLong_cat), function(x) {
  
  # filter the dataset for dietitians with the specific years of practice
  dat_practice <- dat %>% filter(dietitian_forHowLong_cat == x)
  
  results_dat <- data.frame(withinHealthcare_PbdietShouldBeOption  = levels(dat_practice$withinHealthcare_PbdietShouldBeOption),
                            share_yes                   = prop.table(table(dat_practice$withinHealthcare_PbdietShouldBeOption)) %>% as.vector,
                            dietitian_forHowLong_cat    = x) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat_hospitals_years <- plot_dat_list %>% bind_rows() %>% 
  mutate(dietitian_forHowLong_cat    = factor(dietitian_forHowLong_cat, levels = levels(dat$dietitian_forHowLong_cat)),
         withinHealthcare_PbdietShouldBeOption  = factor(withinHealthcare_PbdietShouldBeOption, levels = levels(dat$withinHealthcare_PbdietShouldBeOption)))

# heatmap
plot_dat_hospitals_years %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = withinHealthcare_PbdietShouldBeOption, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["beliefs"]),
                      labels = scales::label_percent()) +
  ggtitle("Should WFPB diets be available in hospitals?\nby dietitians' years of practice") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Barriers research questions

WFPB: Personal barriers

Those that have chosen not to follow a WFPB diet personally, what was the main reasons/barriers in preventing them from adopting one?

Code
plot_dat_personalBarriers <- dat %>% 
  filter(diet_personallyTriedWFPB == "no") %>% 
  select(starts_with("diet_personallyTriedWFPB_no_barriers")) %>% 
  select(-diet_personallyTriedWFPB_no_barriers) %>% 
  pivot_longer(cols = everything()) %>% 
  group_by(name) %>% 
  mutate(sum_yes = sum(value == "yes")) %>% 
  ungroup() %>% 
  arrange(desc(sum_yes)) %>% 
  mutate(name = gsub(pattern = "diet_personallyTriedWFPB_no_barriers", replacement = "", name),
         name = factor(name, levels = rev(unique(name))))
plot_dat_personalBarriers_labels <- plot_dat_personalBarriers %>% 
  group_by(name) %>% 
  summarize(share_yes = first(sum_yes) / n(),
            share_yes_label = paste0(round(100 * share_yes, 0), "%"))

plot_dat_personalBarriers %>% 
  ggplot(aes(y = name)) +
  geom_bar(aes(fill = value), position = "fill") +
  geom_text(data = plot_dat_personalBarriers_labels,
            aes(x = share_yes, label = share_yes_label),
            color = unname(col_vector["barriers"]), nudge_x = 0.07) +
  scale_fill_manual(values = c("lightgray", unname(col_vector["barriers"]))) +
  scale_x_continuous(labels = scales::label_percent()) +
  xlab("Agreement [%]") +
  ggtitle("Barriers - did not personally try WFPB") +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

WFPB: Clients’ barriers

What are some of the barriers their patients have come across in adopting a WFPB diet?

Code
plot_dat_clientsBarriers <- dat %>% 
  select(starts_with("patientsGoWFPB_barriers")) %>% 
  select(-patientsGoWFPB_barriers) %>% 
  pivot_longer(cols = everything()) %>% 
  group_by(name) %>% 
  mutate(sum_yes = sum(value == "yes")) %>% 
  ungroup() %>% 
  arrange(desc(sum_yes)) %>% 
  mutate(name = gsub(pattern = "patientsGoWFPB_barriers", replacement = "", name),
         name = factor(name, levels = rev(unique(name))))
plot_dat_clientsBarriers_labels <- plot_dat_clientsBarriers %>% 
  group_by(name) %>% 
  summarize(share_yes = first(sum_yes) / n(),
            share_yes_label = paste0(round(100 * share_yes, 0), "%"))

plot_dat_clientsBarriers %>% 
  ggplot(aes(y = name)) +
  geom_bar(aes(fill = value), position = "fill") +
  geom_text(data = plot_dat_clientsBarriers_labels,
            aes(x = share_yes, label = share_yes_label),
            color = unname(col_vector["barriers"]), nudge_x = 0.07) +
  scale_fill_manual(values = c("lightgray", unname(col_vector["barriers"]))) +
  scale_x_continuous(labels = scales::label_percent()) +
  xlab("Agreement [%]") +
  ggtitle("Barriers for patients that go WFPB") +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

Do the dietitians think a WFPB diet is financially realisitc for low income clients?

Code
dat %>% 
  ggplot(aes(y = "", fill = WFPB_financiallyRealistic_forLowIncome)) +
  geom_bar(position = "fill") +
  scale_x_continuous("Relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "RdBu") +
  guides(fill = guide_legend(reverse = TRUE)) +
  ggtitle("WFPB diets are financially realistic for low income clients") +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

WFPB: Enough university education?

Did they receive enough education on WFPB diets when in university?

Code
dat %>% 
  ggplot(aes(y = "", fill = education_WFPBeducation)) +
  geom_bar(position = "fill") +
  scale_x_continuous("Relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "RdBu") +
  guides(fill = guide_legend(reverse = TRUE)) +
  ggtitle("Personal WFPB education") +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

WFPB: Prescription confidence

Do they feel confident in prescribing a WFPB diet?

Code
dat %>% 
  ggplot(aes(y = "", fill = WFPBcounselling_confidence)) +
  geom_bar(position = "fill") +
  scale_x_continuous("Relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "RdBu") +
  guides(fill = guide_legend(reverse = TRUE)) +
  ggtitle("Confidence in WFPB counselling") +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

Univariate description

by subspecialty

Code
# prepare the dataset
plot_dat_list <- lapply(1:length(subspecialties_varVector), function(i) {
  
  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # filter the dataset for dietitians with the specific subspecialty
  dat_subsp <- dat
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  results_dat <- data.frame(WFPBcounselling_confidence = levels(dat_subsp$WFPBcounselling_confidence),
             share_yes                   = prop.table(table(dat_subsp$WFPBcounselling_confidence)) %>% as.vector,
             subspecialty                = subsp_label) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat_confidence_sub <- plot_dat_list %>% bind_rows() %>% 
  mutate(subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)),
         WFPBcounselling_confidence = factor(WFPBcounselling_confidence, levels = levels(dat$WFPBcounselling_confidence)))

# heatmap
plot_dat_confidence_sub %>% 
  ggplot(aes(x = subspecialty, y = WFPBcounselling_confidence, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["barriers"]),
                      labels = scales::label_percent()) +
  ggtitle("Confidence in WFPB counselling\nby subspecialties") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

WFPB: Quality of education

Do they feel there is proper education on WFPB diets?

Code
dat %>% 
  ggplot(aes(y = "", fill = education_WFPBeducation)) +
  geom_bar(position = "fill") +
  scale_x_continuous("Relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "RdBu") +
  guides(fill = guide_legend(reverse = TRUE)) +
  ggtitle("Quality of education on WFPB diet during dietetic degree") +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

Univariate description

by subspecialty

Code
# prepare the dataset
plot_dat_list <- lapply(1:length(subspecialties_varVector), function(i) {
  
  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # filter the dataset for dietitians with the specific subspecialty
  dat_subsp <- dat
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  results_dat <- data.frame(education_WFPBeducation = levels(dat_subsp$education_WFPBeducation),
             share_yes                   = prop.table(table(dat_subsp$education_WFPBeducation)) %>% as.vector,
             subspecialty                = subsp_label) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat_eduQuality_sub <- plot_dat_list %>% bind_rows() %>% 
  mutate(subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)),
         education_WFPBeducation = factor(education_WFPBeducation, levels = levels(dat$education_WFPBeducation)))

# heatmap
plot_dat_eduQuality_sub %>% 
  ggplot(aes(x = subspecialty, y = education_WFPBeducation, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["barriers"]),
                      labels = scales::label_percent()) +
  ggtitle("Quality of education on WFPB diet during dietetic degree\nby subspecialties") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

WFPB: Educational resources

Do they feel there are enough educational resources for them to implement a WFPB in their setting?

Code
dat %>% 
  ggplot(aes(y = "", fill = Pbdiet_enoughProperEducationalResources)) +
  geom_bar(position = "fill") +
  scale_x_continuous("Relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "RdBu") +
  guides(fill = guide_legend(reverse = TRUE)) +
  ggtitle("Enough educational resources on WFPB diets?") +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

Univariate description

by subspecialty

Code
# prepare the dataset
plot_dat_list <- lapply(1:length(subspecialties_varVector), function(i) {
  
  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # filter the dataset for dietitians with the specific subspecialty
  dat_subsp <- dat
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  results_dat <- data.frame(Pbdiet_enoughProperEducationalResources = levels(dat_subsp$Pbdiet_enoughProperEducationalResources),
             share_yes                   = prop.table(table(dat_subsp$Pbdiet_enoughProperEducationalResources)) %>% as.vector,
             subspecialty                = subsp_label) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat_eduResources_sub <- plot_dat_list %>% bind_rows() %>% 
  mutate(subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)),
         Pbdiet_enoughProperEducationalResources = factor(Pbdiet_enoughProperEducationalResources, levels = levels(dat$Pbdiet_enoughProperEducationalResources)))

# heatmap
plot_dat_eduResources_sub %>% 
  ggplot(aes(x = subspecialty, y = Pbdiet_enoughProperEducationalResources, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["barriers"]),
                      labels = scales::label_percent()) +
  ggtitle("Enough educational resources on WFPB diets?\nby subspecialties") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

WFPB: Workplace support

Do they feel supported in advocating for a WFPB diet in their workplace?

Code
dat %>% 
  ggplot(aes(y = "", fill = workplace_supportForWFPBadvocacy)) +
  geom_bar(position = "fill") +
  scale_x_continuous("Relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "RdBu") +
  guides(fill = guide_legend(reverse = TRUE)) +
  ggtitle("Workplace support for WFPB advocacy") +
  theme(axis.title.y    = element_blank(),
        legend.title    = element_blank(),
        legend.position = "bottom")

Univariate description

by subspecialty

Code
# prepare the dataset
plot_dat_list <- lapply(1:length(subspecialties_varVector), function(i) {
  
  subsp       <- unname(subspecialties_varVector)[i]
  subsp_label <- names(subspecialties_varVector)[i]
  
  # filter the dataset for dietitians with the specific subspecialty
  dat_subsp <- dat
  colnames(dat_subsp)[colnames(dat_subsp) == subsp] <- "x"
  dat_subsp <- dat_subsp %>% filter(x == "yes")
  
  results_dat <- data.frame(workplace_supportForWFPBadvocacy = levels(dat_subsp$workplace_supportForWFPBadvocacy),
             share_yes                   = prop.table(table(dat_subsp$workplace_supportForWFPBadvocacy)) %>% as.vector,
             subspecialty                = subsp_label) %>% 
    mutate(share_label   = paste0(round(100 * share_yes), "%"))
  
  return(results_dat)
})
plot_dat_workplace_sub <- plot_dat_list %>% bind_rows() %>% 
  mutate(subspecialty = factor(subspecialty, levels = names(subspecialties_varVector)),
         workplace_supportForWFPBadvocacy = factor(workplace_supportForWFPBadvocacy, levels = levels(dat$workplace_supportForWFPBadvocacy)))

# heatmap
plot_dat_workplace_sub %>% 
  ggplot(aes(x = subspecialty, y = workplace_supportForWFPBadvocacy, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("Agreement [%]", low = "white", high = unname(col_vector["barriers"]),
                      labels = scales::label_percent()) +
  ggtitle("Workplace support for WFPB advocacy\nby subspecialties") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank())

Figures for publication

Figure 1: Knowledge score descriptives

Code
fig2_baseSize <- 14

# overall distribution of life cycle knowledge score
meanScore_lifeCycle <- paste0(100 * round(mean(dat$knowledge_score_lifeCycle, na.rm = TRUE), 2), "%")
gg_score_lifeCycle <- ggplot(dat, aes(x = knowledge_score_lifeCycle, y = after_stat(count / sum(count)))) +
  geom_histogram(fill = col_vector["knowledge"], bins = 8) +
  scale_y_continuous("relative frequency", labels = scales::label_percent(),
                     breaks = c(0, .1, .2), minor_breaks = c(.05, .15)) +
  scale_x_continuous("score", labels = scales::label_percent(), limits = c(0,NA),
                     breaks = c(0, .5, 1)) +
  ggtitle("Life cycle", paste("mean score:", meanScore_lifeCycle)) +
  labs(tag = "(A)") +
  theme_minimal(base_size = fig2_baseSize) +
  theme(panel.grid.minor.x = element_blank(),
        panel.grid.major.x = element_blank(),
        plot.tag.position  = c(0,1),
        plot.tag           = element_text(hjust = 0, vjust = 1),
        plot.title         = element_text(hjust = 0.5),
        plot.subtitle      = element_text(hjust = 0.5))

# overall distribution of critical micronutrients knowledge score
meanScore_criticalNutrients <- paste0(100 * round(mean(dat$knowledge_score_criticalNutrients, na.rm = TRUE), 2), "%")
gg_score_criticalNutrients <- ggplot(dat, aes(x = knowledge_score_criticalNutrients, y = after_stat(count / sum(count)))) +
  geom_bar(fill = col_vector["knowledge"]) +
  scale_y_continuous("relative frequency", labels = scales::label_percent(),
                     breaks = c(0, .1, .2), minor_breaks = c(.05, .15, .25)) +
  scale_x_continuous("score", labels = scales::label_percent(),
                     breaks = seq(0, 1, by = .5)) +
  ggtitle("Critical nutrients", paste("mean score:", meanScore_criticalNutrients)) +
  theme_minimal(base_size = fig2_baseSize) +
  theme(panel.grid.minor.x = element_blank(),
        panel.grid.major.x = element_blank(),
        axis.title.y       = element_blank(),
        plot.title         = element_text(hjust = 0.5),
        plot.subtitle      = element_text(hjust = 0.5))

# overall distribution of uncritical micronutrients knowledge score
meanScore_uncriticalNutrients <- paste0(100 * round(mean(dat$knowledge_score_uncriticalNutrients, na.rm = TRUE), 2), "%")
gg_score_uncriticalNutrients <- ggplot(dat, aes(x = knowledge_score_uncriticalNutrients, y = after_stat(count / sum(count)))) +
  geom_bar(fill = col_vector["knowledge"]) +
  scale_y_continuous("relative frequency", labels = scales::label_percent(),
                     breaks = c(0, .2, .4), minor_breaks = c(.1, .3, .5)) +
  scale_x_continuous("score", labels = scales::label_percent(),
                     breaks = seq(0, 1, by = .5), limits = c(0, NA)) +
  ggtitle("Non-critical nutrients", paste("mean score:", meanScore_uncriticalNutrients)) +
  theme_minimal(base_size = fig2_baseSize) +
  theme(panel.grid.minor.x = element_blank(),
        panel.grid.major.x = element_blank(),
        axis.title.y       = element_blank(),
        plot.title         = element_text(hjust = 0.5),
        plot.subtitle      = element_text(hjust = 0.5))

# overall distribution of improved diseases knowledge score
meanScore_diseases <- paste0(100 * round(mean(dat$knowledge_score_diseases, na.rm = TRUE), 2), "%")
gg_score_diseases <- ggplot(dat, aes(x = knowledge_score_diseases, y = after_stat(count / sum(count)))) +
  geom_bar(fill = col_vector["knowledge"]) +
  scale_y_continuous("relative frequency", labels = scales::label_percent(),
                     limits = c(0, .15)) +
  scale_x_continuous("score", labels = scales::label_percent(),
                     breaks = seq(0, 1, by = .5)) +
  ggtitle("Diseases", paste("mean score:", meanScore_diseases)) +
  theme_minimal(base_size = fig2_baseSize) +
  theme(panel.grid.minor   = element_blank(),
        panel.grid.major.x = element_blank(),
        axis.title.y       = element_blank(),
        plot.title         = element_text(hjust = 0.5),
        plot.subtitle      = element_text(hjust = 0.5))


# joint the results tables to create a joint effect plot
results_tab_LCscore <- results_tab_LCscore %>% mutate(model = "life cycle")
results_tab_CNscore <- results_tab_CNscore %>% mutate(model = "critical nutrients")
results_tab_UNscore <- results_tab_UNscore %>% mutate(model = "uncritical nutrients")
results_tab_IDscore <- results_tab_IDscore %>% mutate(model = "diseases")
results_tab <- results_tab_LCscore %>% 
  bind_rows(results_tab_CNscore) %>% 
  bind_rows(results_tab_UNscore) %>% 
  bind_rows(results_tab_IDscore) %>% 
  mutate(model = factor(model, levels = c("life cycle", "critical nutrients",
                                          "uncritical nutrients", "diseases")))

# effect plot
results_tab_edited <- results_tab %>% 
  filter(parameter != "(Intercept)") %>% 
  mutate(parameter = case_when(parameter == "diabetesyes"                  ~ "diabetes",
                               parameter == "elderlyyes"                   ~ "care for the elderly",
                               parameter == "gastroenterologyyes"          ~ "gastroenterology",
                               parameter == "obesityOrWeightManagementyes" ~ "weight management",
                               parameter == "oncologyyes"                  ~ "oncology",
                               parameter == "paediatricsyes"               ~ "paediatrics",
                               parameter == "psychEatingyes"               ~ "eating disorders",
                               grepl("forHowLong", parameter)              ~ "",
                               parameter == "educationPhD"                 ~ "PhD",
                               parameter == "educationPostgraduate degree" ~ "postgraduate degree",
                               parameter == "education_WFPBeducation_num"  ~ "",
                               TRUE                                        ~ parameter),
         parameter = factor(parameter, levels = c("weight management",
                                                  "diabetes",
                                                  "gastroenterology",
                                                  "paediatrics",
                                                  "care for the elderly",
                                                  "oncology",
                                                  "eating disorders",
                                                  "",
                                                  "postgraduate degree",
                                                  "PhD")),
         model = factor(model,
                        levels = c("life cycle", "critical nutrients", "uncritical nutrients", "diseases"),
                        labels = c("life cycle", "critical\nnutrients", "non-critical\nnutrients", "diseases")))
signifEffects_dat <- results_tab_edited %>% 
  filter(pvalue < .05) %>% 
  mutate(y = .175)
gg_effectPlot <- results_tab_edited %>% 
  ggplot(aes(x = parameter)) +
  geom_hline(yintercept = 0, lty = 2, col = "gray50") +
  geom_pointrange(aes(y = estimate, ymin = CI_lower, ymax = CI_upper, col = param_group)) +
  # geom_point(data = signifEffects_dat, aes(y = y), shape = 4, size = 2, stroke = 2,
  #            col = "gray20") +
  geom_point(data = signifEffects_dat, aes(y = y), shape = 8, size = 2, stroke = 0.5,
             col = "gray20") +
  facet_grid(rows = vars(model), cols = vars(param_group), scales = "free_x") +
  scale_y_continuous(limits = c(-.20, .20),
                     breaks = c(-.1, 0, .1),
                     labels = c("-0.1", "0", "+0.1"),
                     minor_breaks = c(-.2, .2)) +
  scale_color_brewer(palette = "Dark2") +
  ggtitle("Knowledge score model results") +
  labs(tag = "(B)") +
  theme_minimal(base_size = fig2_baseSize) +
  theme(axis.text.x        = element_text(angle = 61, hjust = 1, size = 14),
        axis.title.x       = element_blank(),
        strip.text         = element_text(size = 14),
        legend.position    = "none",
        panel.grid.major.x = element_blank(),
        plot.title         = element_text(hjust = 0.5),
        plot.tag.position  = c(0,1),
        plot.tag           = element_text(hjust = 0, vjust = 1),
        plot.margin        = margin(b = 25, t = 15, 5.5, 5.5))


# overall distribution of increased risk of malnutrition / eating disorders
gg_malnutrition <- dat %>% 
  filter(!is.na(Pbdiet_increasedRiskForMalnutrition),
         !is.na(Pbdiet_increasedRiskEatingDisorder)) %>% 
  select(Pbdiet_increasedRiskForMalnutrition, Pbdiet_increasedRiskEatingDisorder) %>% 
  pivot_longer(cols = everything(), names_to = "item") %>% 
  mutate(item  = case_when(item == "Pbdiet_increasedRiskForMalnutrition" ~ "   malnutrition",
                           item == "Pbdiet_increasedRiskEatingDisorder"  ~ "   eating disorders"),
         value = factor(value,
                        levels = c("Strongly disagree", "Disagree", "Not sure", "Agree", "Strongly agree"),
                        labels = c("Strongly disagree", "Disagree", "Not sure", "Agree", "Strongly agree  "))) %>%
  ggplot(aes(y = item, fill = value)) +
  geom_bar(position = "fill") +
  scale_x_continuous("relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "PRGn") +
  guides(fill = guide_legend(reverse = TRUE)) +
  ggtitle("Plant-based diets increase the risk for ...") +
  labs(tag = "(C)") +
  theme_minimal(base_size = fig2_baseSize) +
  theme(axis.title.y      = element_blank(),
        axis.text.y       = element_text(size = 14),
        legend.title      = element_blank(),
        legend.text       = element_text(size = 14, margin = margin(l = 5, r = 15)),
        legend.position   = "bottom",
        panel.grid        = element_blank(),
        plot.title        = element_text(hjust = 0.5),
        plot.tag.position = c(0,1),
        plot.tag          = element_text(hjust = 0, vjust = 1),
        plot.background   = element_rect(fill = "white", color = "white"),
        panel.background  = element_rect(fill = "white", color = "white"))


### joint plot
# 1) align main plots using patchwork
layout <- "
ABCD
EEEE
"
gg_main <- gg_score_lifeCycle + gg_score_criticalNutrients +
  gg_score_uncriticalNutrients + gg_score_diseases + gg_effectPlot +
  patchwork::plot_layout(design = layout, heights = c(.2, .8))

# 2) add the 'increased risk' items using ggpubr since patchwork handles the y axis labels weirdly
ggpubr::ggarrange(plotlist = list(gg_main, gg_malnutrition), nrow = 2, heights = c(.82, .18))

Code
ggsave("../figures/Figure 1 - Knowledge scores.png", width = 10, height = 12.5)

Figure 2: Plant-protein

Code
# overall distribution
gg1 <- dat %>% 
  ggplot(aes(y = "", fill = plantProtein_isIncomplete)) +
  geom_bar(position = "fill") +
  scale_x_continuous("relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "PRGn", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE)) +
  ggtitle("Plant proteins are an incomplete source of protein") +
  labs(tag = "(A)") +
  theme_minimal(base_size  = 12) +
  theme(axis.title.y       = element_blank(),
        legend.title       = element_blank(),
        legend.position    = "bottom",
        legend.text        = element_text(size = 10, margin = margin(l = 5, r = 15)),
        panel.grid.major.y = element_blank(),
        plot.title         = element_text(hjust = 0.5),
        plot.tag.position  = c(0,1),
        plot.tag           = element_text(hjust = 0, vjust = 1))

# heatmap vs. subspecialty
signif_subspecialties <- results_tab_plantProtein %>% 
  filter(pvalue < .05,
         parameter != "(Intercept)") %>% 
  pull(parameter) %>% 
  gsub(pattern = "yes", replacement = "")
gg2 <- plot_dat_proteinSubsp %>% 
  ggplot(aes(x = subspecialty, y = plantProtein_isIncomplete, fill = share_yes)) +
  geom_tile() +
  geom_text(data = plot_dat_proteinSubsp %>%
              filter(!(subspecialty %in% signif_subspecialties)),
            aes(label = share_label), color = "gray40", size = 3) +
  geom_text(data = plot_dat_proteinSubsp %>%
              filter(subspecialty %in% signif_subspecialties),
            aes(label = share_label), color = "gray20", fontface = "bold", size = 3) +
  scale_fill_gradient("agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      labels = scales::label_percent(), limits = c(0, 1)) +
  ggtitle("by areas of specialty") +
  labs(tag = "(B)") +
  theme_minimal(base_size = 12) +
  theme(axis.text.x       = element_text(angle = 45, hjust = 1),
        axis.title        = element_blank(),
        plot.tag.position = c(0,1),
        plot.tag          = element_text(hjust = 0, vjust = 1),
        plot.title         = element_text(hjust = 0.5),
        legend.position   = "none")

# heatmap vs. years of practice
gg3 <- plot_dat_proteinPractice %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = plantProtein_isIncomplete, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      labels = scales::label_percent(), limits = c(0, 1)) +
  ggtitle("by years of practice") +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text.y = element_blank(),
        plot.title         = element_text(hjust = 0.5),
        axis.title  = element_blank())

# joint plot
gg1 / (gg2 + gg3) +
  patchwork::plot_layout(nrow = 2, heights = c(.2,.8))

Code
ggsave("../figures/Figure 2 - Plant protein.png", width = 10, height = 5)

Figure 3: Risk reduction for the three different diseases

Code
update_dietLabels <- function(dat) {
  dat %>% 
    arrange(name) %>% 
    mutate(name = case_when(name == "Medit"       ~ "Mediterranean",
                            name == "National"    ~ "National guidelines",
                            name == "LowCarbs"    ~ "Low carb",
                            name == "LowEnergy"   ~ "Low energy",
                            name == "TDR"         ~ "Total diet replacement",
                            name == "LowFat"      ~ "Low fat",
                            name == "PDR"         ~ "Partial diet replacement",
                            name == "HighProtein" ~ "High protein",
                            name == "Keto"        ~ "Ketogenic",
                            TRUE                  ~ name),
           name = factor(name, levels = unique(as.character(name))))
}
Code
fig3_baseSize <- 14

# overall distribution for type-2 diabetes
plot_dat_labels_T2Diab <- plot_dat_T2Diab %>% 
  filter(name != "NA") %>% 
  group_by(name) %>% 
  summarize(share_yes       = first(sum_yes) / n(),
            share_yes_label = paste0(round(100 * share_yes, 0), "%")) %>% 
  mutate(panel = case_when(name %in% c("WFPB", "Vegan") ~ "Panel A",
                           TRUE                         ~ "Panel B")) %>% 
  update_dietLabels()
ggDiab_overall <- plot_dat_T2Diab %>% 
  update_dietLabels() %>% 
  filter(name != "NA") %>% 
  mutate(panel = case_when(name %in% c("WFPB", "Vegan") ~ "Panel A",
                           TRUE                         ~ "Panel B")) %>% 
  ggplot(aes(y = name)) +
  geom_bar(aes(fill = value), position = "fill") +
  geom_text(data = plot_dat_labels_T2Diab, aes(x = share_yes, label = share_yes_label),
            nudge_x = 0.09, color = unname(col_vector["knowledge"])) +
  scale_fill_manual("diet\nrecommendation", values = c("lightgray", unname(col_vector["knowledge"]))) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_x_continuous("agreement [%]", labels = scales::label_percent()) +
  facet_grid(rows = vars(panel), scales = "free_y", space = "free_y") +
  ggtitle("Type-2 diabetes") +
  labs(tag = "(A)") +
  theme_minimal(base_size = fig3_baseSize) +
  theme(axis.title.y      = element_blank(),
        axis.text.x       = element_blank(),
        axis.title.x      = element_blank(),
        axis.text.y       = element_text(size = fig3_baseSize),
        panel.grid        = element_blank(),
        legend.text       = element_text(size = fig3_baseSize),
        strip.text        = element_blank(),
        plot.tag.position = c(0,1),
        plot.tag          = element_text(hjust = 0, vjust = 1),
        plot.title        = element_text(hjust = 0.5))


# overall distribution for cardiovascular disease
plot_dat_labels_cardiovascDisease <- plot_dat_cardiovascDisease %>% 
  filter(name != "NA") %>% 
  group_by(name) %>% 
  summarize(share_yes       = first(sum_yes) / n(),
            share_yes_label = paste0(round(100 * share_yes, 0), "%")) %>% 
  mutate(panel = case_when(name %in% c("WFPB", "Vegan") ~ "Panel A",
                           TRUE                         ~ "Panel B")) %>% 
  update_dietLabels()
ggCardio_overall <- plot_dat_cardiovascDisease %>% 
  update_dietLabels() %>% 
  filter(name != "NA") %>% 
  mutate(panel = case_when(name %in% c("WFPB", "Vegan") ~ "Panel A",
                           TRUE                         ~ "Panel B")) %>% 
  ggplot(aes(y = name)) +
  geom_bar(aes(fill = value), position = "fill") +
  geom_text(data = plot_dat_labels_cardiovascDisease, aes(x = share_yes, label = share_yes_label),
            nudge_x = 0.09, color = unname(col_vector["knowledge"])) +
  scale_fill_manual("diet\nrecommendation", values = c("lightgray", unname(col_vector["knowledge"]))) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_x_continuous("agreement [%]", labels = scales::label_percent()) +
  facet_grid(rows = vars(panel), scales = "free_y", space = "free_y") +
  ggtitle("Cardiovascular disease") +
  labs(tag = "(B)") +
  theme_minimal(base_size = fig3_baseSize) +
  theme(axis.title.y      = element_blank(),
        axis.text.x       = element_blank(),
        axis.text.y       = element_text(size = fig3_baseSize),
        axis.title.x      = element_blank(),
        legend.text       = element_text(size = fig3_baseSize),
        panel.grid        = element_blank(),
        strip.text        = element_blank(),
        plot.tag.position = c(0,1),
        plot.tag          = element_text(hjust = 0, vjust = 1),
        plot.title        = element_text(hjust = 0.5))

# overall distribution for weight loss
plot_dat_labels_weightLoss <- plot_dat_weightLoss %>% 
  filter(name != "NA") %>% 
  group_by(name) %>% 
  summarize(share_yes       = first(sum_yes) / n(),
            share_yes_label = paste0(round(100 * share_yes, 0), "%")) %>% 
  mutate(panel = case_when(name %in% c("WFPB", "Vegan") ~ "Panel A",
                           TRUE                         ~ "Panel B")) %>% 
  update_dietLabels()
ggWeight_overall <- plot_dat_weightLoss %>% 
  update_dietLabels() %>% 
  filter(name != "NA") %>% 
  mutate(panel = case_when(name %in% c("WFPB", "Vegan") ~ "Panel A",
                           TRUE                         ~ "Panel B")) %>% 
  ggplot(aes(y = name)) +
  geom_bar(aes(fill = value), position = "fill") +
  geom_text(data = plot_dat_labels_weightLoss, aes(x = share_yes, label = share_yes_label),
            nudge_x = 0.09, color = unname(col_vector["knowledge"])) +
  scale_fill_manual("diet\nrecommendation", values = c("lightgray", unname(col_vector["knowledge"]))) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_x_continuous("agreement [%]", labels = scales::label_percent()) +
  facet_grid(rows = vars(panel), scales = "free_y", space = "free_y") +
  ggtitle("Weight loss") +
  labs(tag = "(C)") +
  theme_minimal(base_size = fig3_baseSize) +
  theme(axis.title.y      = element_blank(),
        axis.text.x       = element_text(size = fig3_baseSize),
        axis.text.y       = element_text(size = fig3_baseSize),
        panel.grid        = element_blank(),
        legend.text       = element_text(size = fig3_baseSize),
        strip.text        = element_blank(),
        plot.tag.position = c(0,1),
        plot.tag          = element_text(hjust = 0, vjust = 1),
        plot.title        = element_text(hjust = 0.5))

# use identical axis limits for the heatmaps
axis_limits <- c(0,1)

# subspecialty-specific distribution for type-2 diabetes
# model_T2D_WFPB$results_tab %>% 
#   filter(pvalue < .05,
#          parameter != "(Intercept)")
ggDiab_sub <- plot_dat_T2Diab_sub %>% 
  filter(diet != "NA") %>% 
  mutate(panel = case_when(diet %in% c("WFPB", "Vegan") ~ "Panel A",
                           TRUE                         ~ "Panel B")) %>% 
  ggplot(aes(x = subspecialty, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  facet_grid(rows = vars(panel), scales = "free_y", space = "free_y") +
  ggtitle("by areas of specialty") +
  theme_minimal(base_size = fig3_baseSize) +
  theme(axis.text  = element_blank(),
        axis.title = element_blank(),
        panel.grid = element_blank(),
        strip.text = element_blank(),
        plot.title = element_text(hjust = 0.5))

# subspecialty-specific distribution for cardiovascular disease
# model_CD_WFPB$results_tab %>%
#   filter(pvalue < .05,
#          parameter != "(Intercept)")
ggCardio_sub <- plot_dat_cardiovascDisease_sub %>% 
  filter(diet != "NA") %>% 
  mutate(panel = case_when(diet %in% c("WFPB", "Vegan") ~ "Panel A",
                           TRUE                         ~ "Panel B")) %>% 
  ggplot(aes(x = subspecialty, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  facet_grid(rows = vars(panel), scales = "free_y", space = "free_y") +
  theme_minimal(base_size = fig3_baseSize) +
  theme(axis.title = element_blank(),
        axis.text  = element_blank(),
        panel.grid = element_blank(),
        strip.text = element_blank())

# subspecialty-specific distribution for weight loss
signif_subspecialties <- model_WL_WFPB$results_tab %>%
  filter(pvalue < .05,
         parameter != "(Intercept)") %>% 
  mutate(parameter = case_when(grepl("WeightManag", parameter) ~ "weight management",
                               TRUE                            ~ parameter)) %>% 
  pull(parameter) %>% 
  gsub(pattern = "yes", replacement = "")
plot_dat_weightLoss_sub_final <- plot_dat_weightLoss_sub %>% 
  filter(diet != "NA") %>% 
  mutate(panel = case_when(diet %in% c("WFPB", "Vegan") ~ "Panel A",
                           TRUE                         ~ "Panel B"))
ggWeight_sub <- ggplot(plot_dat_weightLoss_sub_final,
                       aes(x = subspecialty, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(data = plot_dat_weightLoss_sub_final %>% 
              filter((diet != "WFPB") | !(subspecialty %in% signif_subspecialties)),
            aes(label = share_label), color = "gray40", size = 3) +
  geom_text(data = plot_dat_weightLoss_sub_final %>% 
              filter((diet == "WFPB") & (subspecialty %in% signif_subspecialties)),
            aes(label = share_label), color = "gray20", fontface = "bold", size = 3) +
  scale_fill_gradient("agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  facet_grid(rows = vars(panel), scales = "free_y", space = "free_y") +
  theme_minimal(base_size = fig3_baseSize) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = fig3_baseSize),
        axis.text.y = element_blank(),
        axis.title  = element_blank(),
        panel.grid  = element_blank(),
        strip.text  = element_blank())

# 'years of practice'-specific distribution for type-2 diabetes
ggDiab_lon <- plot_dat_T2Diab_lon %>% 
  filter(diet != "NA") %>% 
  mutate(panel = case_when(diet %in% c("WFPB", "Vegan") ~ "Panel A",
                           TRUE                         ~ "Panel B")) %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  facet_grid(rows = vars(panel), scales = "free_y", space = "free_y") +
  ggtitle("by years of practice") +
  theme_minimal(base_size = fig3_baseSize) +
  theme(axis.title = element_blank(),
        axis.text  = element_blank(),
        panel.grid = element_blank(),
        strip.text = element_blank(),
        plot.title = element_text(hjust = 0.5))

# 'years of practice'-specific distribution for cardiovascular disease
ggCardio_lon <- plot_dat_cardiovascDisease_lon %>% 
  filter(diet != "NA") %>% 
  mutate(panel = case_when(diet %in% c("WFPB", "Vegan") ~ "Panel A",
                           TRUE                         ~ "Panel B")) %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  facet_grid(rows = vars(panel), scales = "free_y", space = "free_y") +
  theme_minimal(base_size = fig3_baseSize) +
  theme(axis.title = element_blank(),
        axis.text  = element_blank(),
        panel.grid = element_blank(),
        strip.text = element_blank())

# 'years of practice'-specific distribution for weight loss
ggWeight_lon <- plot_dat_weightLoss_lon %>% 
  filter(diet != "NA") %>% 
  mutate(panel = case_when(diet %in% c("WFPB", "Vegan") ~ "Panel A",
                           TRUE                         ~ "Panel B")) %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = diet, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("agreement [%]", low = "white", high = unname(col_vector["knowledge"]),
                      limits = axis_limits, labels = scales::label_percent()) +
  facet_grid(rows = vars(panel), scales = "free_y", space = "free_y") +
  theme_minimal(base_size = fig3_baseSize) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = fig3_baseSize),
        axis.text.y = element_blank(),
        axis.title  = element_blank(),
        panel.grid  = element_blank(),
        strip.text  = element_blank())

# joint plot
(ggDiab_overall + ggDiab_sub + ggDiab_lon) /
  (ggCardio_overall + ggCardio_sub + ggCardio_lon) /
  (ggWeight_overall + ggWeight_sub + ggWeight_lon) +
  patchwork::plot_layout(guides = "collect")

Code
ggsave("../figures/Figure 3 - Risk reduction.png", width = 12, height = 13)

Figure 4: Beliefs

Code
fig4_baseSize <- 14

# overall distribution of long-term sustainability
ggSus_overall <- dat %>% 
  ggplot(aes(y = "", fill = WFPB_isSustainableLongTerm)) +
  geom_bar(position = "fill") +
  scale_x_continuous("relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer("agreement", type = "div", palette = "BrBG") +
  guides(fill = guide_legend(reverse = TRUE)) +
  ggtitle("Long-term sustainability of WFPB diets") +
  labs(tag = "(A)") +
  theme_minimal(base_size = fig4_baseSize) +
  theme(axis.title.y      = element_blank(),
        panel.grid        = element_blank(),
        plot.title        = element_text(hjust = 0.5),
        plot.tag.position = c(0,1),
        plot.tag          = element_text(hjust = 0, vjust = 1))

# subspecialty-specific distribution of long-term sustainability
signif_subspecialties <- results_tab_longTermSus %>%
  filter(pvalue < .05,
         parameter != "(Intercept)",
         !grepl("education", parameter)) %>% 
  pull(parameter) %>% 
  gsub(pattern = "yes", replacement = "")
ggSus_sub <- plot_dat_sustainability_sub %>% 
  ggplot(aes(x = subspecialty, y = WFPB_isSustainableLongTerm, fill = share_yes)) +
  geom_tile() +
  geom_text(data = plot_dat_sustainability_sub %>% 
              filter(!(subspecialty %in% signif_subspecialties)),
            aes(label = share_label), color = "gray40", size = 3) +
  geom_text(data = plot_dat_sustainability_sub %>% 
              filter(subspecialty %in% signif_subspecialties),
            aes(label = share_label), color = "gray20", fontface = "bold", size = 3) +
  scale_fill_gradient("agreement [%]", low = "white", high = unname(col_vector["beliefs"]),
                      labels = scales::label_percent(), limits = c(0,1)) +
  ggtitle("by areas of specialty") +
  theme_minimal(base_size = fig4_baseSize) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank(),
        plot.title  = element_text(hjust = 0.5),
        panel.grid  = element_blank())

# 'years of practice'-specific distribution of long-term sustainability
ggSus_years <- plot_dat_sustainability_years %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = WFPB_isSustainableLongTerm, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("agreement [%]", low = "white", high = unname(col_vector["beliefs"]),
                      labels = scales::label_percent(), limits = c(0,1)) +
  ggtitle("by years of practice") +
  theme_minimal(base_size = fig4_baseSize) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text.y = element_blank(),
        axis.title  = element_blank(),
        plot.title  = element_text(hjust = 0.5),
        panel.grid  = element_blank())

# overall distribution of availability in hospitals
ggHosp_overall <- dat %>% 
  ggplot(aes(y = "", fill = withinHealthcare_PbdietShouldBeOption)) +
  geom_bar(position = "fill") +
  scale_x_continuous("relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer("Agreement", type = "div", palette = "BrBG") +
  ggtitle("      Should WFPB diets be available in hospitals?") +
  labs(tag = "(B)") +
  theme_minimal(base_size = fig4_baseSize) +
  theme(axis.title.y      = element_blank(),
        legend.position   = "none",
        panel.grid        = element_blank(),
        plot.title        = element_text(hjust = 0.5),
        plot.tag.position = c(0,1),
        plot.tag          = element_text(hjust = 0, vjust = 1))

# subspecialty-specific distribution of availability in hospitals
ggHosp_sub <- plot_dat_hospitals_sub %>% 
  ggplot(aes(x = subspecialty, y = withinHealthcare_PbdietShouldBeOption, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("agreement [%]", low = "white", high = unname(col_vector["beliefs"]),
                      labels = scales::label_percent(), limits = c(0,1)) +
  ggtitle("by areas of specialty") +
  theme_minimal(base_size = fig4_baseSize) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text.y = element_blank(),
        axis.title  = element_blank(),
        plot.title  = element_text(hjust = 0.5),
        panel.grid  = element_blank())

# 'years of practice'-specific distribution of availability in hospitals
ggHosp_years <- plot_dat_hospitals_years %>% 
  ggplot(aes(x = dietitian_forHowLong_cat, y = withinHealthcare_PbdietShouldBeOption, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 3) +
  scale_fill_gradient("agreement [%]", low = "white", high = unname(col_vector["beliefs"]),
                      labels = scales::label_percent(), limits = c(0,1)) +
  ggtitle("by years of practice") +
  theme_minimal(base_size = fig4_baseSize) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text.y = element_blank(),
        axis.title  = element_blank(),
        plot.title  = element_text(hjust = 0.5),
        panel.grid  = element_blank())

# joint plot
layout <- "
AABB
CDEF
"
ggSus_overall + ggHosp_overall +
  ggSus_sub + ggSus_years + ggHosp_sub + ggHosp_years +
  patchwork::plot_layout(nrow = 2, heights = c(.2, .8), design = layout, guides = "collect")

Code
ggsave("../figures/Figure 4 - Beliefs.png", width = 15, height = 5)

Figure 5: Barriers

Code
fig5_baseSize <- 18

# overall distribution of personal barriers
plot_dat_personalBarriers_labels <- plot_dat_personalBarriers_labels %>% 
  mutate(name = case_when(name == "Dairy"          ~ "Excluding dairy",
                          name == "Eggs"           ~ "Excluding eggs",
                          name == "Micronutrients" ~ "Micronutrient\nconcerns",
                          name == "Protein"        ~ "Protein inadequacy",
                          name == "PBdiet"         ~ "Perceived difficulty\nto follow PB diet",
                          name == "Cultural"       ~ "Adapting ethnic\nor culture foods",
                          name == "Healthy"        ~ "Not healthy",
                          TRUE                     ~ name),
         name = factor(name, levels = name))
gg_personalBarriers_overall <- plot_dat_personalBarriers %>% 
  mutate(name = case_when(name == "Dairy"          ~ "Excluding dairy",
                          name == "Eggs"           ~ "Excluding eggs",
                          name == "Micronutrients" ~ "Micronutrient\nconcerns",
                          name == "Protein"        ~ "Protein inadequacy",
                          name == "PBdiet"         ~ "Perceived difficulty\nto follow PB diet",
                          name == "Cultural"       ~ "Adapting ethnic\nor culture foods",
                          name == "Healthy"        ~ "Not healthy",
                          TRUE                     ~ name)) %>%
  ggplot(aes(y = name)) +
  geom_bar(aes(fill = value), position = "fill") +
  geom_text(data = plot_dat_personalBarriers_labels,
            aes(x = share_yes, label = share_yes_label),
            color = unname(col_vector["barriers"]), nudge_x = 0.1, size = 6) +
  scale_fill_manual(values = c("lightgray", unname(col_vector["barriers"]))) +
  scale_x_continuous("agreement [%]", labels = scales::label_percent()) +
  ggtitle("Barriers for personal change                         ") +
  labs(tag = "(A)") +
  theme_minimal(base_size = fig5_baseSize) +
  theme(axis.title.y      = element_blank(),
        legend.title      = element_blank(),
        legend.position   = "none",
        panel.grid.major.y = element_blank(),
        panel.grid.minor   = element_blank(),
        plot.tag.position = c(0,1),
        plot.tag          = element_text(hjust = 0, vjust = 1),
        plot.title        = element_text(hjust = 0.5))

# overall distribution of clients' barriers
plot_dat_clientsBarriers_labels <- plot_dat_clientsBarriers_labels %>% 
  mutate(name = case_when(name == "Difficulty"  ~ "Perceived difficulty\nto follow a PB diet",
                          name == "MealPrep"    ~ "Meal preparation",
                          name == "Interest"    ~ "Lack of interest",
                          name == "FoodGroups"  ~ "Excluding\nfood groups",
                          name == "Knowledge"   ~ "Lack of\nhealth knowledge",
                          name == "Beliefs"     ~ "Negative beliefs",
                          name == "Ethnic"      ~ "Adapting ethnic\nor culture foods",
                          name == "Protein"     ~ "Protein inadequacy",
                          name == "FoodOptions" ~ "Food options",
                          TRUE                     ~ name),
         name = factor(name, levels = name))
gg_clientsBarriers_overall <- plot_dat_clientsBarriers %>% 
  mutate(name = case_when(name == "Difficulty"  ~ "Perceived difficulty\nto follow a PB diet",
                          name == "MealPrep"    ~ "Meal preparation",
                          name == "Interest"    ~ "Lack of interest",
                          name == "FoodGroups"  ~ "Excluding\nfood groups",
                          name == "Knowledge"   ~ "Lack of\nhealth knowledge",
                          name == "Beliefs"     ~ "Negative beliefs",
                          name == "Ethnic"      ~ "Adapting ethnic\nor culture foods",
                          name == "Protein"     ~ "Protein inadequacy",
                          name == "FoodOptions" ~ "Food options",
                          TRUE                     ~ name)) %>%
  ggplot(aes(y = name)) +
  geom_bar(aes(fill = value), position = "fill") +
  geom_text(data = plot_dat_clientsBarriers_labels,
            aes(x = share_yes, label = share_yes_label),
            color = unname(col_vector["barriers"]), nudge_x = 0.1, size = 6) +
  scale_fill_manual(values = c("lightgray", unname(col_vector["barriers"]))) +
  scale_x_continuous("agreement [%]", labels = scales::label_percent()) +
  ggtitle("Barriers for patients' change                         ") +
  theme_minimal(base_size  = fig5_baseSize) +
  theme(axis.title.y       = element_blank(),
        legend.title       = element_blank(),
        legend.position    = "none",
        panel.grid.major.y = element_blank(),
        panel.grid.minor   = element_blank(),
        plot.tag.position  = c(0,1),
        plot.tag           = element_text(hjust = 0, vjust = 1),
        plot.title         = element_text(hjust = 0.5))

# use identical legend limits for subspecialty-specific plots
axis_limits_barriers <- c(0, 1)

# overall distribution of the confidence prescribing WFPB diets
gg_confidence_overall <- dat %>% 
  filter(WFPBcounselling_confidence != "Not applicable") %>% 
  mutate(WFPBcounselling_confidence = factor(WFPBcounselling_confidence,
                                             levels = c("Not confident",
                                                        "Slightly confident",
                                                        "Somewhat confident",
                                                        "Confident",
                                                        "Fairly confident",
                                                        "Completely confident"),
                                             labels = c("Not confident",
                                                        "Slightly confident",
                                                        "Somewhat confident    ",
                                                        "Confident",
                                                        "Fairly confident",
                                                        "Completely confident    "))) %>%
  ggplot(aes(y = "", fill = WFPBcounselling_confidence)) +
  geom_bar(position = "fill") +
  scale_x_continuous("relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "RdBu") +
  guides(fill = guide_legend(nrow = 3, reverse = TRUE)) +
  ggtitle("    Confidence in WFPB counselling") +
  labs(tag = "(C)") +
  theme_minimal(base_size = fig5_baseSize) +
  theme(axis.title.y      = element_blank(),
        legend.title      = element_blank(),
        legend.position   = "bottom",
        legend.text       = element_text(margin = margin(l = 5, r = 15)),
        plot.tag.position = c(-.05,1),
        plot.tag          = element_text(hjust = 0, vjust = 1),
        plot.title        = element_text(hjust = 0.5))
gg_confidence_legend <- ggpubr::get_legend(gg_confidence_overall) %>% 
  ggpubr::as_ggplot()
gg_confidence_overall <- gg_confidence_overall +
  theme(legend.position = "none")

# subspecialty-specific distribution of the confidence prescribing WFPB diets
gg_confidence_sub <- plot_dat_confidence_sub %>% 
  ggplot(aes(x = subspecialty, y = WFPBcounselling_confidence, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 4) +
  scale_fill_gradient("agree-\nment\n[%]", low = "white", high = unname(col_vector["barriers"]),
                      limits = axis_limits_barriers, labels = scales::label_percent()) +
  ggtitle("by areas of specialty") +
  theme_minimal(base_size = fig5_baseSize) +
  theme(axis.text.x = element_blank(),
        axis.title  = element_blank(),
        plot.title  = element_text(hjust = 0.5),
        panel.grid  = element_blank())

# overall distribution of education quality on WFPB
gg_eduQuality_overall <- dat %>% 
  ggplot(aes(y = "", fill = education_WFPBeducation)) +
  geom_bar(position = "fill") +
  scale_x_continuous("relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "RdBu") +
  guides(fill = guide_legend(nrow = 2, reverse = TRUE)) +
  ggtitle("        Quality of education on WFPB diet") +
  labs(tag = "(D)") +
  theme_minimal(base_size = fig5_baseSize) +
  theme(axis.title.y      = element_blank(),
        legend.title      = element_blank(),
        legend.position   = "bottom",
        legend.text       = element_text(margin = margin(l = 5, r = 15)),
        plot.tag.position = c(-.05,1),
        plot.tag          = element_text(hjust = 0, vjust = 1),
        plot.title        = element_text(hjust = 0.5))
gg_eduQuality_legend <- ggpubr::get_legend(gg_eduQuality_overall) %>% 
  ggpubr::as_ggplot()
gg_eduQuality_overall <- gg_eduQuality_overall +
  theme(legend.position = "none")

# subspecialty-specific distribution of education quality on WFPB
gg_eduQuality_sub <- plot_dat_eduQuality_sub %>% 
  ggplot(aes(x = subspecialty, y = education_WFPBeducation, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 4) +
  scale_fill_gradient("agree-\nment\n[%]", low = "white", high = unname(col_vector["barriers"]),
                      limits = axis_limits_barriers, labels = scales::label_percent()) +
  ggtitle("by areas of specialty") +
  theme_minimal(base_size = fig5_baseSize) +
  theme(axis.text.x = element_blank(),
        axis.title  = element_blank(),
        plot.title  = element_text(hjust = 0.5),
        panel.grid  = element_blank())

# overall distribution of educational resources on WFPB
gg_eduResources_overall <- dat %>% 
  ggplot(aes(y = "", fill = Pbdiet_enoughProperEducationalResources)) +
  geom_bar(position = "fill") +
  scale_x_continuous("relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "RdBu") +
  guides(fill = guide_legend(nrow = 2, reverse = TRUE)) +
  ggtitle("                        Enough educational resources on WFPB diets") +
  labs(tag = "(E)") +
  theme_minimal(base_size = fig5_baseSize) +
  theme(axis.title.y      = element_blank(),
        legend.title      = element_blank(),
        legend.position   = "bottom",
        legend.text       = element_text(margin = margin(l = 5, r = 15)),
        plot.tag.position = c(-.05,1),
        plot.tag          = element_text(hjust = 0, vjust = 1),
        plot.title        = element_text(hjust = 0.5))
gg_eduResources_legend <- ggpubr::get_legend(gg_eduResources_overall) %>% 
  ggpubr::as_ggplot()
gg_eduResources_overall <- gg_eduResources_overall +
  theme(legend.position = "none")

# subspecialty-specific distribution of educational resources on WFPB
gg_eduResources_sub <- plot_dat_eduResources_sub %>% 
  ggplot(aes(x = subspecialty, y = Pbdiet_enoughProperEducationalResources, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 4) +
  scale_fill_gradient("agree-\nment\n[%]", low = "white", high = unname(col_vector["barriers"]),
                      limits = axis_limits_barriers, labels = scales::label_percent()) +
  ggtitle("by areas of specialty") +
  theme_minimal(base_size = fig5_baseSize) +
  theme(axis.text.x = element_blank(),
        axis.title  = element_blank(),
        plot.title  = element_text(hjust = 0.5),
        panel.grid  = element_blank())

# overall distribution of workplace support
gg_workplace_overall <- dat %>% 
  mutate(workplace_supportForWFPBadvocacy = factor(workplace_supportForWFPBadvocacy,
                                                   levels = c("Strongly unsupported",
                                                              "Somewhat unsupported",
                                                              "Neutral",
                                                              "Somewhat supported",
                                                              "Strongly supported"),
                                                   labels = c("Strongly unsupported",
                                                              "Somewhat unsupported",
                                                              "Neutral",
                                                              "Somewhat supported    ",
                                                              "Strongly supported"))) %>% 
  ggplot(aes(y = "", fill = workplace_supportForWFPBadvocacy)) +
  geom_bar(position = "fill") +
  scale_x_continuous("relative frequency [%]", labels = scales::label_percent()) +
  scale_fill_brewer(type = "div", palette = "RdBu") +
  guides(fill = guide_legend(nrow = 3, reverse = TRUE)) +
  ggtitle("              Workplace support for WFPB advocacy") +
  labs(tag = "(F)") +
  theme_minimal(base_size = fig5_baseSize) +
  theme(axis.title.y      = element_blank(),
        legend.title      = element_blank(),
        legend.position   = "bottom",
        legend.margin     = margin(l = -20),
        legend.text       = element_text(margin = margin(l = 5, r = 10)),
        plot.tag.position = c(-.05,1),
        plot.tag          = element_text(hjust = 0, vjust = 1),
        plot.title        = element_text(hjust = 0.5))
gg_workplace_legend <- ggpubr::get_legend(gg_workplace_overall) %>% 
  ggpubr::as_ggplot()
gg_workplace_overall <- gg_workplace_overall +
  theme(legend.position = "none")

# subspecialty-specific distribution of workplace support
gg_workplace_sub <- plot_dat_workplace_sub %>% 
  ggplot(aes(x = subspecialty, y = workplace_supportForWFPBadvocacy, fill = share_yes)) +
  geom_tile() +
  geom_text(aes(label = share_label), color = "gray40", size = 4) +
  scale_fill_gradient("agree-\nment\n[%]", low = "white", high = unname(col_vector["barriers"]),
                      limits = axis_limits_barriers, labels = scales::label_percent()) +
  ggtitle("by areas of specialty") +
  theme_minimal(base_size = fig5_baseSize) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title  = element_blank(),
        plot.title  = element_text(hjust = 0.5),
        panel.grid  = element_blank())


# joint plot
layout <- "
ACK
ADK
AEL
AFL
BGM
BHM
BIN
BJN
"
gg_personalBarriers_overall + gg_clientsBarriers_overall +
  gg_confidence_overall + gg_confidence_legend +
  gg_eduQuality_overall + gg_eduQuality_legend + 
  gg_eduResources_overall + gg_eduResources_legend + 
  gg_workplace_overall + gg_workplace_legend +
  gg_confidence_sub + gg_eduQuality_sub + gg_eduResources_sub + gg_workplace_sub +
  patchwork::plot_layout(design = layout, heights = c(3,7,3,7,3,7,3,7), widths = c(.3,.4,.3),
                         guides = "collect")

Code
ggsave("../figures/Figure 5 - Barriers.png", width = 19, height = 14)

Supplementary Figure 3: Disease management

Code
sFig2_baseSize <- 14

dat %>% 
  select(starts_with("Pbdiet_improvedConditions_")) %>% 
  pivot_longer(cols = everything()) %>% 
  group_by(name) %>% 
  mutate(sum_yes = sum(value == "yes")) %>% 
  ungroup() %>% 
  arrange(desc(sum_yes)) %>% 
  mutate(name = gsub(pattern = "Pbdiet_improvedConditions_", replacement = "", name),
         name = case_when(name == "heartDisease"     ~ "Heart disease",
                          name == "highCholesterol"  ~ "High cholesterol",
                          name == "hypertension"     ~ "Hypertension",
                          name == "T2Diabetes"       ~ "T2DM",
                          name == "obesity"          ~ "Obesity",
                          name == "certainCancers"   ~ "Certain cancers",
                          name == "stroke"           ~ "Stroke",
                          name == "fattyLiver"       ~ "Fatty liver disease",
                          name == "depression"       ~ "Depression",
                          name == "alzheimer"        ~ "Alzheimer's",
                          name == "vascularDementia" ~ "Vascular dementia",
                          name == "chronicKidney"    ~ "Chronic kidney disease",
                          TRUE                       ~ name),
         name = factor(name, levels = rev(unique(name)))) %>% 
  ggplot(aes(y = name, fill = value)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("lightgray", unname(col_vector["knowledge"]))) +
  scale_x_continuous("Agreement [%]", labels = scales::label_percent()) +
  ylab("Chronic disease") +
  ggtitle("Improved conditions through plant-based diet") +
  theme_minimal(base_size = sFig2_baseSize) +
  theme(legend.title     = element_blank(),
        legend.position  = "bottom",
        plot.title       = element_text(hjust = 0.5),
        plot.background  = element_rect(fill = "white", color = "white"),
        panel.background = element_rect(fill = "white", color = "white"),
        panel.grid       = element_blank())

Code
ggsave("../figures/Supplementary Figure 3 - Disease management.png", width = 8, height = 5)